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Abstract 


A  systems  approach  by  which  to  understand  the  characteristics  and  scaling  of  wall  turbulence,  resolvent 
analysis  (or  critical  later  framework)  has  been  developed,  implemented  and  reported  in  a  range  of  archival 
publications,  beginning  with  McKeon  and  Sharma  [2010].  The  approach  is  in  essence  a  transfer  function 
analysis  of  the  linear  Navier-Stokes  operator  in  the  streamwise/spanwise  wavenumber  and  temporal  frequency 
space,  which  identifies  and  makes  use  of  the  low  rank  nature  of  this  operator  in  order  to  achieve  order 
reduction.  The  outputs  have  been  shown  to  capture  key  elements  of  the  turbulence,  from  statistics  to  spectral 
information  to  arrays  of  coherent  structure.  Supporting  experiments  and  modeling  have  been  performed  to 
provide  first  comprehensive  evidence  of  where  in  wavenumber-frequency  space  turbulence  activity  occurs  by 
characterizing  the  three-dimensional  spectrum  of  the  streamwise  velocity  fluctuations  using  time-resolved 
particle  image  velocimetry  at  a  range  of  wall-normal  locations,  and  extending  a  two-dimensional,  three- 
component  model  to  provide  an  approximate  input  mean  velocity  profile  to  the  analysis.  The  full  resolvent 
analysis  has  the  potential  to  provide  a  self-consistent  and  self-sustaining  model  of  turbulence,  in  which 
the  velocity  outputs  have  the  correct  amplitude  and  phase  to  support  the  appropriate  input  forcing  and 
(presently  assumed)  mean  velocity  profile.  Further  work  is  proposed  which  exploits  sparsity  in  the  frequency 
domain  associated  with  the  resolvent  formulation  to  further  reduce  the  complexity  of  the  representation  of 
turbulence  and  to  continue  to  explore  the  implications  for  the  scaling  of  wall  turbulence. 
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1  Introduction 


The  research  performed  under  this  PECASE  award  targeted  multi-scale  experiments  and  modeling  in  wall 
turbulence.  The  major  contribution  of  the  work  has  been  the  introduction  and  development  of  a  novel, 
systems  approach  to  understanding  the  behavior  of  unperturbed  wall  turbulence,  which  we  refer  to  as 
“resolvent  analysis” . 


1.1  Objectives 

The  specific  objectives  of  the  research,  as  listed  in  the  original  proposal,  were  as  follows. 

1.  Perform  first  experimental  measurements  of  the  full  frequency-streamwise/spanwise  wavenumber  spec¬ 
trum  in  canonical  and  non-canonical  turbulent  boundary  layers  using  fast  Particle  Image  Velocimetry 
in  order  to  identify  the  influence  of  multi-scale  structure,  wall-normal  location  and  Reynolds  number 
on  the  inferred  dispersion  relationship  relating  frequency  and  wavenumber; 

2.  Develop  novel  non-normal  mode  analysis  of  harmonic  forcing  of  a  canonical  wall-bounded  turbulent 
flow  in  order  to  categorize  the  extent,  location  and  scaling  properties  of  structures  identified  using  a 
singular  value  decomposition  of  the  linear  operator; 

3.  Concatenate  the  results  of  the  experimental  and  analytical  projects  in  order  to  take  steps  towards  a 
novel,  low  order,  yet  multi-scale,  model  of  wall  turbulence. 

Significant  effort  was  spent  on  each  of  these  objectives,  and  we  believe  that  progress  has  been  made  in  each. 

The  results  have  been  disseminated  widely.  In  this  report,  the  over-arching  approach  will  be  summarized  in 

the  following  section,  before  further  technical  details  will  be  given  in  each  area. 

1.2  Dissemination  of  this  work  (excluding  numerous  oral  presentations) 

1.2.1  Archival  journals 

1.  McKeon,  B.  J.  and  Sharma,  A.  (2010)  “A  critical  layer  framework  for  turbulent  pipe  flow”,  J.  Fluid 
Mech.,658,  336-382. 

2.  Smits,  A.  J.,  McKeon,  B.  J.  and  Marusic,  I.  “High  Reynolds  number  wall  turbulence”,  Annual  Review 
of  Fluid  Mechanics,  43,  2011. 

3.  LeHew,  J.,  Guala,  M.  and  McKeon,  B.  J.  “A  study  of  the  three-dimensional  spectral  energy  distribution 
in  a  zero  pressure  gradient  turbulent  boundary  layer”,  Expts.  in  Fluids,  10.1007/s00348-011-1117-z 
(2011). 

4.  Jacobi,  I.  and  McKeon,  B.  J.  “New  perspectives  on  the  impulsive  roughness  perturbation  of  a  turbulent 
boundary  layer”,  J.  Fluid  Mech.,  677,  179-203  (2011). 

5.  Bourguignon,  J.-L.  and  McKeon,  B.  J.  “A  streamwise-constant  model  of  turbulent  pipe  flow”,  Phys. 
Fluids,  23,  095111  (2011). 

6.  Sharma,  A.  and  McKeon,  B.  J.  “On  coherent  structure  in  wall  turbulence”,  J.  Fluid  Mech.  728,  196-238 
(2013);  also  ArXiV1301.7580,  2013. 

7.  LeHew,  J.,  Guala,  M.  and  McKeon,  B.  J.  “Time-resolved  measurements  of  coherent  structures  in  the 
turbulent  boundary  layer”,  Expts.  in  Fluids,  54(4)  (2013). 

8.  McKeon,  B.  J.,  Jacobi,  I  and  Sharma,  A.  S.  “Experimental  manipulation  of  wall  turbulence:  a  systems 
approach”,  Phys.  Fluids  25(3)  10.1063/1.4793444  (2013). 
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9.  Jacobi,  I.  and  McKeon,  B.  J.  “Phase  relationships  between  large  and  small  scales  in  the  turbulent 
boundary  layer”,  Expts.  in  Fluids,  10.1007/s00348-013-1481-y  (2013). 

10.  Moarref,  R.,  Sharrna,  A.  S.,  Tropp,  J.  A.  and  McKeon,  B.  J.  “Model-based  scaling  of  the  streamwise 
energy  density  in  high-Reynolds-number  turbulent  channels”,  J.  Fluid  Mech.,  734,  275-316  (2013). 

11.  Moarref,  R.,  Jovanovic,  M.  R.,  Tropp,  J.  A.,  Slrarma,  A.  S.  and  McKeon,  B.  J.  “A  low-order  decom¬ 
position  of  turbulent  channel  flow  via  resolvent  analysis  and  convex  optimization”,  Phys.  Fluids,  26, 
051701  (2014);  also  ArXiV  1401.6417,  2014. 

12.  Gomez  Carrasco,  F.,  Blackburn,  H.,  Rudrnan,  M.,  McKeon,  B.,  Luhar,  M.,  Moarref,  R.  and  Sharrna, 
A.  “On  the  origin  of  frequency  sparsity  in  direct  numerical  simulations  of  turbulent  pipe  flow” ,  Phys. 
Fluids  26,  101703  (2014). 

13.  Moarref,  R.,  Sharrna,  A.  S.,  Tropp,  J.  A.  and  McKeon,  B.  J.  “A  foundation  for  analytical  developments 
in  the  logarithmic  region  of  turbulent  channels”  (under  review). 

1.2.2  Conference  papers 

1.  LeHew,  J.,  Guala,  M.  and  McKeon,  B.  J.  “A  study  of  convective  velocities  in  a  zero  pressure  gradient 
turbulent  boundary  layer”  AIAA-2010-4474. 

2.  McKeon,  B.  J.  and  Sharrna,  A.  S.,  “Energetic  modes  in  turbulent  pipe  flow  from  resolvent  analysis” 
AIAA-2010-697. 

3.  LeHew,  J.,  Guala,  M.  and  McKeon,  B.  J.  “Structural  and  spectral  analyisis  of  the  turbulent  boundary 
layer  using  time-resolved  2DPIV”  PIV-11,  Kobe,  Japan  (Ninth  International  Symposium  on  Particle 
Image  Velocimetry),  July  2011. 

4.  Bourguignon,  J.-L.  and  McKeon,  B.  J.  “A  streamwise-constant  model  for  turbulent  pipe  flow”  TSFP-7, 
Ottowa  (Turbulent  Shear  Flow  Phenomena  Symposium  -7),  July  2011. 

5.  LeHew,  J.,  Guala,  M.  and  McKeon,  B.  J.  “A  simultaneous  physical  and  spectral  analyisis  of  the 
turbulent  boundary  layer  via  time-resolved  particle  image  velocimetry”  TSFP-7,  Ottowa  (Turbulent 
Shear  Flow  Phenomena  Symposium  -7),  July  2011. 

6.  Jacobi,  I.  and  McKeon,  B.  J.  “Critical-layer  behavior  in  an  impulsively-perturbed  turbulent  boundary 
layer”  TSFP-7,  Ottowa  (Turbulent  Shear  Flow  Phenomena  Symposium  -7),  July  2011. 

7.  Slrarma,  A.  and  McKeon,  B.  J.  “Velocity  statistics  and  structure  in  pipe  turbulence  derived  from  a 
simple  critical-layer  model”  TSFP-7,  Ottowa  (Turbulent  Shear  Flow  Phenomena  Symposium  -7),  July 
2011. 

8.  Slrarma,  A.  S.  and  McKeon,  B.  J.  “Closing  the  loop:  an  explicit  treatment  of  the  nonlinearity  in  the 
resolvent  analysis  of  wall  turbulence”  AIAA-2013-3118. 

9.  Moarref,  R.,  Sharrna,  A.,  Tropp,  J.  and  McKeon,  B.  J.  “Prediction  of  the  streamwise  turbulent  energy 
spectrum  in  high  Reynolds  number  channel  flows”  AIAA-2013-2480. 

10.  Moarref,  R.,  Sharrna,  A.  S.,  Tropp,  J.  A.  and  McKeon,  B.  J.  “A  predictive  low-rank  model  for  turbu¬ 
lent  energy  spectra  at  high  Reynolds  numbers”  TSFP-8,  Poitiers  (Turbulent  Shear  Flow  Phenomena 
Symposium  -8),  August  2013. 

11.  Moarref,  R.,  Sharrna,  A.,  Tropp,  J.  and  McKeon,  B.  J.  “On  scaling  and  interaction  of  the  geometrically 
self-similar  modes  in  turbulent  channels”  USNCTAM  2014. 

12.  Moarref,  R.,  Park,  J.-S.,  Sharnra,  A.  S.,  Willis,  A.  P.,  Graham,  M.  D.  and  McKeon,  B.  J.  “Approx¬ 
imation  of  the  exact  traveling  wave  solutions  in  wall-bounded  flows  using  resolvent  modes”  TSFP-9, 
Melbourne  (Turbulent  Shear  Flow  Phenomena  Symposium  -9),  July  2015  (to  appear). 
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2  Approach  and  Summary  of  Results 

2.1  Set-up 

Our  approach  to  the  multiscale  nature  of  wall  turbulence  has  been  to  develop  a  scale-by-scale  analysis 
with  its  foundations  in  systems  analysis  and  connections  with  the  non-normal  mechanisms  of  disturbance 
amplification  first  explored  in  the  extensive  literature  transient  growth.  This  section  gives  a  high-level 

description  of  the  approach;  further  technical  details  on  each  aspect  can  be  found  in  the  subsequent  sections 

and  the  publications  listed  in  section  1. 

While  many  studies  have  identified  the  potential  for  large  amplification  of  disturbances  due  to  the  non¬ 
normality  of  the  linearized  Navier-Stokes  operator  (LNS),  these  have  mostly  been  framed  in  the  context 
of  the  search  for  the  globally  most  amplified  disturbances  and  associated  flow  response.  In  addition,  the 
majority  of  work  has  been  in  the  time  domain,  considering  transient  growth.  These  studies  uniformly 

return  streamwise  constant,  kx  -  0,  disturbances  as  the  worst-case,  most  amplified  inputs  and  outputs.  The 

connection  to  observable,  “real”  multiscale  turbulence  then  requires  further  modeling,  which  has  had  limited 
success  to-date.  Ad-hoc  assumptions  are  required  to  “shape”  input  disturbances  to  approach  outputs  that 
have  physical  relevance.  The  relevant  literature  is  reviewed  in  section  3. 

In  the  present  work,  we  have  sought  to  reframe  the  problem  in  terms  of  the  disturbances  giving  rise  to 
realizable  features  of  turbulence,  with  a  view  to  reconciling  a  range  of  observations  and  hypotheses  concerning 
its  structure  and  origin.  As  such,  the  analysis  is  constructed  in  the  space  in  which  experimental  and  numerical 
results  are  often  reported,  namely  a  Fourier  decomposition  of  the  velocity  field.  In  particular,  we  consider  the 
triple  Fourier  decomposition  in  the  wall-parallel  directions,  and  in  time,  leading  to  disturbance  amplification 
systems  for  each  assembly  of  streamwise  and  spanwise  wavenumbers,  and  temporal  frequency.  The  picture 
at  each  wavenumber/frequency  combination  is  essentially  that  sketched  in  figure  1,  in  which  the  lower  block 
represents  the  linear  dynamics  of  fluctuations  around  the  turbulent  mean  velocity  profile,  forced  by  nonlinear 
triadic  interactions  from  other  wavenumbers  and  frequencies. 

At  the  heart  of  the  theory  developed  under  this  grant  are  two  insights.  Firstly,  that  turbulence  is  robust. 
Secondly,  that  the  robustness  comes  from  a  feedback  loop  involving  the  Navier-Stokes  nonlinearity,  which 
is  energy-conserving,  and  a  highly  selective  linear  system,  which  leads  to  highly  directional  amplification. 
That  is,  the  nonlinearity  provides  the  correct  shaping  of  the  input  forcing  to  sustain  an  output  of  a  velocity 
field  that  corresponds  to  real  turbulence. 

Robustness  may  be  defined  as  the  ability  of  a  system  to  operate  in  a  similar  way  under  a  wide  range 
of  conditions  and  uncertainties.  Wall  turbulence  is  nothing  if  not  robust,  indeed,  it  is  the  persistence  of 
complex,  but  recognizable,  features  that  characterizes  these  flows  as  turbulent.  Essentially,  then,  in  seeking 
to  understand  the  processes  that  characterize  wall  turbulence,  we  are  seeking  to  explain  its  robustness.  This  is 
not  to  say  that  “turbulence  is  linear” ,  but  rather  that  the  processes  that  drive  turbulence  and  select  structure 
do  appear  to  have  a  strong  linear  flavor  (as  hinted  at  by  previous  studies) .  We  deliberately  confine  nonlinear 
effects  are  confined  to  the  determination  of  the  mean  velocity  profile  and  a  conservative  nonlinearity.  By  thus 
doing,  we  are  able  to  make  significant  progress  by  consideration  of  the  linear  (not  linearized)  Navier-Stokes 
operator.  We  approach  the  turbulent  flow  as  a  system,  to  give  us  information  about  how  the  system  behaves 
as  a  whole,  inclusive  of  the  nonlinearity  and  with  a  view  to  making  direct  connections  to  experimental  and 
numerical  observations  of  wall  turbulence.  A  further  benefit  of  the  systems  paradigm  advocated  here  is  that 
closed-loop  control  thinking  and  technology  will  carry  over  in  the  future  for  control  designs,  offering  the 
opportunity  for  sophisticated  closed- loop  strategies  (effectively  based  on  a  perfect-system  model),  different 
from  the  open-loop  “manipulation”  that  has  been  common  in  the  literature. 

The  full  nonlinear  behavior  of  a  turbulent  flow  is  probably  too  complicated  to  model  in  a  simple  manner, 
because  the  interconnection  between  scales  governed  by  the  nonlinear  term  would  require  modeling  all  scales 
at  once.  One  approach,  then,  is  to  pretend  we  know  nothing  of  the  nonlinear  interactions  between  scales 
and  concentrate  on  the  part  we  can  most  easily  analyze,  the  linearized  system.  However,  in  our  case,  the 
nonlinearity  is  not  assumed  away;  instead,  we  treat  it  as  an  unstructured  forcing  that  has  already  acted  to 
support  an  (assumed)  turbulent  mean  profile  of  the  appropriate  form.  The  resulting  decomposition  of  the 
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Figure  1:  A  high-level  description  of  the  turbulence  process.  The  lower  block  contains  the  linear  dynamics 
of  the  fluctuations  interacting  with  the  mean  velocity  profile. 

NSE  as  a  linear  system  driven  by  an  unknown  nonlinear  term  works  because  the  linearized  system  is  so 
selective  that  the  exact  form  of  the  forcing  is  almost  unimportant  under  this  formulation. 

In  keeping  with  similar  approaches  as  applied  to  other  flows  in  the  literature,  we  refer  to  our  analysis  as 
resolvent  analysis.  At  the  schematic  level,  this  analysis,  reported  in  detail  in  McKeon  and  Sharma  [2010] 
and  section  4,  can  be  represented  as  in  figure  2.  It  provides,  first,  an  exact  representation  of  the  NSE  in 
which  the  output  from  a  linear  sub-system  at  each  wavenumber /frequency  combination  interacts  with  all 
others  to  provide  forcing  to  all  linear  sub-systems.  The  operator  in  each  linear  sub-system,  the  resolvent, 
describes  the  transfer  of  energy  into  the  turbulent  fluctuations,  along  with  one  describing  the  sustenance 
of  the  mean  profile  that  appears  in  all  the  resolvents.  The  inter-wavenumber  energy  transfer  then  follows 
from  the  nonlinear  term  as  a  consequence  of  gradients  in  energy  distribution  and  triadic  spectral  interaction 
mechanisms.  In  this  sense,  the  resolvent  formulation  serves  as  a  road  map  for  the  energy  transfer  in  wall 
turbulence,  and  captures  key  dynamical  processes. 

Such  information  makes  possible  theoretical  extensions  to  the  basic  framework,  in  order  to  predict  the 
influence  of  prescribed  external  inputs  and  to  design  forcing  strategies  for  specific,  application-optimized 
modification  of  turbulent  spectra.  Secondly,  the  formulation  admits  significant  mathematical  simplification 
associated  with  the  low- rank  nature  of  the  resolvent,  which  can  be  exploited  to  enable  conceptually  and 
numerically  simple  interrogation,  replication  and  extension  of  known  results  in  wall  turbulence,  ranging 
from  the  statistical  to  structural. 

Figure  2  shows  the  full  schematic  of  the  resolvent  analysis,  including  the  singular  value  decomposition  of 
the  resolvent,  which  constitutes  a  gain-based  decomposition  in  the  wall-normal  direction.  The  first  singular 
input  (forcing)  and  output  (  response)  modes  have  been  the  focus  of  our  investigation  because  they  capture 
many  aspects  of  wall  turbulence.  The  figure  identifies  other  aspects  of  the  analysis  requiring  additional 
study,  namely  the  mechanisms  for  generating  the  mean  velocity  profile,  U(y ),  which  is  assumed  from  exper¬ 
imental  data  in  the  present  analysis,  and  the  required  characteristics  of  the  nonlinear  term  itself,  which  can 
be  expressed  in  terms  of  determining  the  appropriate  response  mode  amplitudes  and  phases  contributing 
to  nonlinear  interactions.  Both  aspects  require  explicit  treatment  of  nonlinearity,  and  resolution  of  both  is 
required  for  the  system  to  be  self-sustaining,  i.e.  to  fully  capture  the  dynamics  of  wall  turbulence.  Pre¬ 
liminary  steps  towards  both  of  these  outstanding  issues  have  been  taken  under  this  grant.  Time-resolved 
particle  image  velocimetry  (TRPIV)  has  been  used  to  provide  a  first  quantification  of  the  three-dimensional 
(streamwise/spanwise  wavenumber  and  temporal  frequency)  spectrum  of  the  streamwise  velocity  through¬ 
out  the  turbulent  boundary  layer  and  to  identify  convection  velocities  associated  with  individual  coherent 
vortical  turbulent  structures,  and  the  so-called  “2D/3C”  (two-dimensional,  three  velocity  component)  model 
for  the  mean  velocity  profile  has  been  extended  from  plane  Couette  flow  to  pipe  (Hagen-Poiseuille)  flow. 

Note  that  the  ability  to  interact  with  the  resolvent  model  via  external  forcing  of  one  transfer  function, 
i.e.  actuation  at  a  single  wavenumber/frequency  combination,  was  introduced  under  AFOSR  grant  9550- 
08-1-0049  “Turbulent  boundary  layer  response  to  perturbation  by  dynamically  rough  surfaces”  and  is  under 
continued  exploration  through  a  Basic  Research  Initiative  grant  9550-12-1-0479  “BRI  Wall  turbulence  with 
designer  properties:  identification,  characterization  &  manipulation  of  dominant  energy  pathways” .  The 
latter  grant  has  also  investigated  the  pressure  signal  associated  with  individual  resolvent  modes  and  means 
by  which  to  implement  linear  control  of  the  resolvent  system  using  wall  actuation. 

The  detailed  outcomes  from  the  present  grant  are  summarized  in  sections  4-6.  A  summary  of  the  key  results 
is  given  below. 


Figure  2:  Resolvent  analysis  as  an  interconnected  system  of  transfer  functions,  'H{kx,kz,L>j).  Each  transfer 
function  is  decomposed  into  the  most  “dangerous”,  or  most  amplified  directions  using  a  singular  value 
decomposition,  such  that  input  forcing,  /  is  ranked  at  each  wavenumber/frequency  triplet  based  on  gain 
(singular,  value  a)  and  giving  rise  to  output  u. 


2.2  Summary  of  results 

2.2.1  Key  features  of  wall  turbulence  captured  by  resolvent  analysis 

Our  publications  [McKeon  and  Sharma,  2010,  Sharma  and  McKeon,  2013b,  Moarref  et  al.,  2013]have  de¬ 
veloped  several  aspects  of  wall  turbulence  that  are  captured,  extended  and/or  explained  using  resolvent 
analysis.  The  key  observations  are  summarized  here. 

Examination  of  the  form  of  the  first  singular  response  modes  (henceforth  velocity  response  modes)  reveals 
that  these  in  isolation  replicate  the  streak/roll  form  that  is  commonly  observed  in  wall  turbulence.  Figure  3 
shows  the  velocity  structure  associated  with  a  wavenumber/frequency  combination  representative  of  the 
near-wall  turbulence  cycle,  capturing  the  wall-normal  extent  known  to  be  associated  with  this  turbulent 
activity.  The  critical  layer  behavior  revealed  by  the  resolvent  analysis  ensures  that  this  structure  is  observed 
at  all  scales,  both  attached  to  and  detached  from  the  wall  (in  the  sense  of  whether  the  velocity  footprint 
reaches  the  wall  or  not). 

Examination  of  the  swirling  (vortical)  motion  associated  with  individual  response  modes  led  to  the  identifi¬ 
cation  of  a  periodic  array  of  prograde  and  retrograde  hairpin  vortices  (with  heads  of  both  senses  of  rotation, 
where  prograde  indicates  a  sense  aligned  with  the  mean  spanwise  vorticity  and  retrograde  indicates  the 
opposite  sense)  associated  with  each  mode  attached  to  the  wall.  Whereas  the  velocity  response  modes  can 
be  linearly  summed,  the  swirling  diagnostic  is  nonlinear,  and  the  mean  velocity  was  shown  to  preferentially 
emphasize  prograde  vortices  while  suppressing  retorgrade  ones,  in  line  with  experimental  observations.  Fur¬ 
ther,  linear  addition  of  even  three  response  modes  at  different  wavenumber/frequency  combinations  (chosen 
to  be  triadically  consistent,  such  that  there  is  potential  for  the  modes  to  self-excite  through  the  nonlinear 
forcing  term),  leads  to  a  complicated  assembly  of  hairpin  vortices,  which  we  deemed  “an  ideal  packet”  with 
characteristics  consistent  with  observed  hairpin  packets,  see  figure  4.  The  resolvent  analysis  allows  a  rela¬ 
tively  simple  origin  of  the  complex  “forests  of  hairpins”  observed  in  high  Reynolds  number  simulations  of 
wall  turbulence,  as  a  linear  assembly  of  velocity  response  modes  giving  rise  to  a  nonlinear  swirling  field. 

In  subsequent  work  Moarref  et  al.  [2013],  geometrical  self-similarity  of  the  resolvent  in  the  logarithmic  region 
has  been  identified,  giving  rise  to  self-similar  hierarchies  of  vortices  as  shown  in  figure  5.  A  paper  on  further 
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Figure  3:  Shape  of  the  first  singular  mode  representative  of  the  dominant  near  wall  motions.  Color 
denotes  isosurfaces  of  streamwise  velocity  (streaks) ,  where  red  and  blue  correspond  to  high  and  low  velocity 
respectively  relative  to  the  mean  flow  (heading  into  the  page),  and  the  white  arrows  show  the  sense  of  the 
in-plane  velocity  field. 
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Figure  4:  Isosurfaces  of  constant  swirling  strength  (33%  of  maximum  value)  for  the  “ideal  packet”  at 
R+  -  1800,  color-coded  by  the  local  (model)  azimuthal  vorticity.  Red  and  blue  denote  rotation  in  and 
counter  to  the  sense  of  the  classical  hairpin  vortex  or  prograde  and  retrograde  vortices,  respectively. 
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Figure  5:  Eddy  structure  at  two  scales  within  a  self-similar  hierarchy  identified  using  resolvent  analysis. 


developments  associated  with  this  self-similarity,  specifically  with  regards  to  self-similar  interaction  between 
hierarchies  forcing  other  hierarchies  is  currently  under  review,  [Moarref  et  ah,  2014b]. 

In  Moarref  et  al.  [2013],  it  was  shown  that  self-similarity  of  the  mean  velocity  profile,  could  be  used  to 
identify  regions  of  self-similarity  of  the  resolvent  and  forcing  magnitudes,  such  that  the  variance  of  the 
velocity  fluctuations  could  be  modeled  satisfactorily  to  high  (atmospheric)  Reynolds  numbers.  An  example 
of  the  comparison  between  experimental/numerical  results  and  the  resolvent  model  for  channel  flow  is  shown 
in  figure  6.  In  further  work  Moarref  et  al.  [2014a],  we  have  shown  that  the  full  velocity  spectra  can  be  well- 
matched  (figure  7)  by  formulation  of  the  appropriate  optimization  problem  and  using  up  to  twelve  resolvent 
modes  at  each  wavenumber/frequency  combination  active  in  turbulence.  This  constitutes  a  model  with  less 
than  0.03%  of  the  degrees  of  freedom  of  the  full  simulation  on  which  the  resolvent  modes  were  projected. 
The  agreement  is  excellent  for  the  streamwise  velocity,  which  dominates  the  energy  norm  used  in  the  singular 
value  decomposition  of  the  resolvent,  but  work  continues  on  the  relatively  weaker  wall-normal  velocity. 

This  work  also  supported  the  earlier  observation  [Bourguignon  et  ah,  2014b]  of  sparsity  in  the  frequency 
domain  when  DNS  data  was  subjected  to  the  three-dimensional  ( k,n,uj )  Fourier  transform.  In  a  recent 
collaborative  publication  [Gomez  Carrasco  et  al.,  2014]  we  have  shown  that  this  sparsity  is  a  real  and 
significant  feature  of  numerical  simulations  that  is  directly  attributable  to  the  spatial  discretization  employed. 
Figure  8  shows  a  comparison  of  the  peak  gain  (amplification,  singular  value)  over  a  range  of  azimuthal 
wavenumbers  and  temporal  frequencies  for  the  continuous  case  and  in  the  presence  of  spatial  discretization. 
The  “peaky”  structure  of  the  discretized  case,  which  is  reflected  as  a  sparsity  in  the  frequency  domain, 
suggests  a  further  order  reduction  for  future  models. 
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Figure  6:  Comparison  of  experimental  and  resolvent  model  results  for  the  streamwise  energy  intensity. 
Red,  blue,  green  lines  show  model  results  at  ReT  =  3165,4000,5813;  symbols  represent  experimental  results 
for  ReT  =  3165,4000  (❖,  v)  [Monty,  2005]  and  ReT  -  4000,6000  (o,  x)  [Schultz  and  Flack,  2013].  Note  the 
influence  of  spatial  averaging  on  the  experimental  profiles  close  to  the  wall,  i.e.  the  disagreement  between 
model  and  experiment  appears  to  be  dominated  by  experimental  bias. 


Figure  7:  Time-averaged  streamwise  velocity  spectra:  from  the  DNS  of  Hoyas  and  Jimenez  [2006]  for 
ReT  =  2003  (solid  lines),  model-based  spectra  with  the  optimal  weights  using  N  =  2,6,12  (a-c)  resolvent 
modes  per  wavenumber/frequency  triplet  and  0  <  c  <  Uci.  Contours  show  10%  to  90%  of  the  maximum  in 
the  DNS  data  with  increments  of  20%. 
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Figure  8:  Distribution  of  resolvent  amplification  (singular  value)  in  the  energetically  active  subset  of  az¬ 
imuthal  wavenumbers  n  and  frequencies  u>/ 2tt  for  the  first  singular  vector  at  Rer  =  314.  (a)  One-dimensional 
resolvent  model  assuming  continuum  of  wavenumbers  available;  (b)  two-dimensional  resolvent  model  incor¬ 
porating  spatial  discretization  adopted  in  DNS. 
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Figure  9:  Premultiplied  3D  spectrum  of  the  streamwise  velocity  fluctuations  in  a  turbulent  boundary  layer. 
Equal  area  within  an  isocontour  corresponds  to  equal  energy. 


2.2.2  TRPIV  characterization  of  the  three-dimensional  power  spectrum  of  streamwise  veloc¬ 
ity  fluctuations 

The  above  description  of  resolvent  analysis  underscores  the  importance  of  understanding  the  distribution  of 
turbulent  energy  in  the  three-dimensional  wavenumber/frequency  space,  in  order  to  provide  targeted  input. 
It  is  an  interesting  problem  in  its  own  right,  since  both  experimental  and  numerical  data  are  presented  after 
integration  in  one  or  two  of  those  dimensions:  typical  experimental  diagnostics  either  implicitly  average  in 
one  dimension  (e.g.  hot  wires)  or  are  used  in  a  mode  that  averages  over  one  dimension  in  order  to  reduce  the 
amount  of  data  storage  required.  One  notable  exception  is  the  work  of  Morrison  and  Kronauer  [1969]  and 
Morrison  et  al.  [1971],  who  obtained  three-dimensional  information  in  a  wall- normal  plane  by  an  exhaustive 
study  that  involved  changing  the  separation  between  two  hot  wires  in  order  to  obtain  spatial  information.  A 
typical  means  to  move  between  wavenumber  and  frequency  space  is  Taylor’s  hypothesis  of  frozen  turbulence, 
in  which  the  local  mean  velocity  is  used  as  the  convection  velocity,  such  that  U(y)  =  to/k,  but  there  are 
known  issues  with  this  assumption,  particularly  for  the  large  scale  motions  close  to  the  wall.  The  response 
mode  shapes  illuminated  by  the  resolvent  analysis  underscore  the  reason  for  these  discrepancies,  namely  that 
the  critical  layers  for  the  larger,  attached  scales  lie  far  from  the  wall,  while  their  footprints  reach  down  close 
to  the  wall  where  the  local  velocity  is  much  less  than  the  correct  convection  velocity  at  the  critical  layer. 

Under  this  grant,  and  published  in  two  papers  in  Experiments  in  Fluids,  LeHew  et  al.  [2011]  and  LeHew  et  al. 
[2013],  direct  measurements  of  the  fully  time-  and  space- resolved  velocity  field  in  wall-parallel  planes  were 
obtained  using  a  5kHz,  time-resolved  PIV  (TRPIV)  system  in  a  turbulent  boundary  layer  developed  over  a 
flat  plate  in  a  water  channel.  This  information  was  processed  to  obtain  the  full  three-dimensional  spectrum 
of  the  streamwise  velocity  fluctuations  in  regions  corresponding,  approximately,  to  the  buffer,  logarithmic 
and  wake  regions  of  the  boundary  layer.  An  example  of  a  premultiplied  three-dimensional  spectrum  obtained 
near  the  top  of  the  logarithmic  region  is  given  in  figure  9.  The  elongated  shape  in  the  to  -  k  plane  indicates 
that  Taylor’s  hypothesis  is  a  good  first  order  estimate  this  far  from  the  wall;  any  asymmetries  in  the  deviation 
from  a  straight  line  correspond  to  failures  of  the  hypothesis  in  an  average  sense. 
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(a)  Histogram  comparison  (b)  Cartoon  of  inferred  structures 

Figure  10:  (a)  Comparison  of  convection  velocity  histograms  at  y/6  =  0.08  and  y/6  =  0.30  leading  to  the 
identification  of  3  distinct  regions.  The  labeled  regions  are  referenced  in  (b)  showing  the  structure  shapes 
inferred  from  the  histograms  in  (a). 


2.2.3  TRPIV  tracking  and  characterization  of  individual  turbulent  structures  and  their  con¬ 
vection  velocities 

The  time-resolved  PIV  experiments  above  were  also  used  to  study  coherent  structures  in  wall-parallel  planes, 
to  track  them  and  identify  structure-specific  trajectory  and  convection  velocity  statistics,  and  to  propose  a 
model  for  the  types  of  structures  observed  by  linking  observations  at  different  wall-normal  heights.  Figure  10a 
shows  the  convection  velocity  “signatures”  of  “swirling  coherent  structures”  (SCS),  which  are  identified 
by  identifying  swirling  motion  in  the  wall-parallel  plane  and  are  considered  to  be  a  proxy  for  coherent 
vortical  motions  with  greater  wall-normal  extent.  The  left-most  curve  corresponds  to  the  distribution  of 
convection  velocities  observed  in  the  buffer  region  measurement  plane  and  the  (overlapping)  rightmost  curve 
to  the  velocities  in  the  wake  region.  Statistically,  velocities  that  are  populated  in  both  planes  are  likely  to 
correspond  to  tall  structures  with  a  signature  in  both  planes,  as  shown  in  schematic  in  figure  lOab.  The 
detailed  observations  in  this  work  were  shown  to  be  consistent  with  Townsend’s  attached  eddy  hypothesis, 
as  identified  directly  from  experimental  data. 


2.2.4  Model  for  the  mean  velocity 

In  the  resolvent  analysis  reported  here,  a  mean  velocity  profile  is  required  as  input  to  the  linear  resol¬ 
vent  operator.  Experimental  data  from  the  Princeton  Superpipe  McKeon  et  al.  [2004]  was  used  for  pipe 
flow[McKeon  and  Sharma,  2010],  and  an  eddy  viscosity  formulation  for  channel  flow  [Moarref  et  al.,  2013]. 
Another  model  for  the  mean  velocity  was  explored  under  this  grant,  namely  an  extension  of  earlier  work 
concerning  the  development  of  the  turbulent  mean  profile  in  Couette  flow  [Gayme  et  al.,  2010]  to  pipe  flow, 
Bourguignon  and  McKeon  [2011].  The  spirit  of  this  investigation  was  to  provide  a  rough  starting  condition 
for  the  resolvent  analysis,  knowing  that  the  correct  assembly  of  mean  profile  and  resolvent  modes  would  be 
self-sustaining,  i.e.  it  would  return  the  correct  mode  amplitudes  and  phases  and  mean  velocity. 

The  so-called  “2D/3C”  model  uses  either  a  model  streamfunction  in  the  cross-stream  plane  or  stochastic 
noise  as  input  to  the  equation  governing  the  streamwise-averaged  velocity  profile,  and  continues  to  return 
encouraging  agreement  with  real  mean  velocity  profiles.  Figure  11  shows  the  mean  velocity  profiles  (b) 
arising  from  the  (v,w)  forcing  streamfunctions  shown  in  (a).  The  agreement  with  DNS  is  promising,  even 
for  such  simplified  streamfunctions,  indicating  a  degree  of  robustness  of  the  mean  velocity  profile  in  addition 
to  the  turbulent  velocities.  In  Bourguignon  and  McKeon  [2011],  stochastic  forcing  was  also  used  in  place  of 
the  simple  streamfunctions  and  led  to  replication  of  several  phenomena  observed  during  pipe  flow  transition. 
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Figure  11:  (a)Streamfunctions  ^ri>0-c(»7)  and  (b)  corresponding  velocity  profiles  Uo(r])  for  = 

0.033(77  -  3r]3  +  2t74)  (thin  solid),  H/i  5(77)  =  0.7(7?  -  3r/3  +  2t ;4)2  (dashed),  and  =  14(77  _  37y3  +  27?4)3 

(dash-dot)  and  experimental  velocity  profile  of  den  Toonder  and  Nieuwstadt  [1997]  at  Re  =  24600  (thick 
solid). 


Preliminary  studies  of  the  resolvent  using  approximate  mean  velocity  profiles  indicated  a  relative  insensitivity 
to  deviations  from  the  real  profile,  with  the  exception  of  regions  of  high  shear,  i.e.  close  to  the  wall.  As  such, 
approximation  of  the  mean  velocity  profile  by  a  model  such  as  the  2D/3C  approach  are  likely  to  constitute 
a  reasonable  starting  condition  for  a  fully-closed  resolvent  analysis.  Such  a  study,  hopefully  leading  to  a 
self-sustaining  model  of  wall  turbulence  in  due  course,  is  reserved  for  future  work,  however  the  results  under 
this  grant  provide  a  promising  first  step. 
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3  Background 


A  full  review  of  the  state  of  the  art  in  wall  turbulence  was  co-authored  in  2011  Snrits  et  al.  [2011]  under 
this  grant.  A  brief  background  and  literature  review  relevant  to  the  work  described  later  in  the  report  is 
given  below.  Further  background  information  can  be  obtained  by  reference  to  the  associated  archival  journal 
papers. 

3.1  A  brief  review  of  linear  processes  in  wall  turbulence 

In  recent  years,  dynamical  systems  approaches  have  had  intriguing  success  in  predicting  the  features  of 
transition,  particularly  in  (linearly  stable)  pipe  flow  [Kerswell,  2005],  as  reviewed  by  Eckhardt  et  al.  [2007]. 
Our  understanding  of  transition  has  also  benefited  from  a  model  in  which  small  disturbances  develop  over 
infinite  or  finite  time  horizons  [Trefethen  et  al.,  1993,  Butler  and  Farrell,  1992a,  Schmid  and  Henningson, 
2001].  In  the  latter  picture,  the  system  is  explored  by  investigating  the  properties  of  the  linearized  Navier- 
Stokes  system,  via  either  its  linear  stability  (eigenvalues)  or  more  recently  by  a  transient  growth  analysis 
(which  investigates  perturbation  growth  on  a  finite  instead  of  infinite  time  horizon).  These  analyses  are 
linear  in  nature,  though  often  they  go  on  to  consider  secondary  instabilities  arising  from  the  nonlinearity. 

One  problem  with  applying  perturbation  analysis  to  turbulence  is  the  selection  of  an  appropriate  base  flow 
around  which  to  linearize,  since  the  turbulent  mean  velocity  profile  does  not  constitute  a  solution  of  the 
Navier-Stokes  equations  (NSE)  and  the  transport  due  to  turbulent  fluctuations  must  be  treated  explicitly. 
However  there  have  been  notable  efforts  to  extend  this  analysis  to  characterize  stability  and  transient  growth 
in  turbulence,  including  the  work  of  Butler  and  Farrell  [1992b],  del  Alamo  and  Jimenez  [2006],  Cossu  et  al. 
[2009]  and  Hwang  and  Cossu  [2010],  who  used  a  turbulent  eddy  viscosity  to  circumvent  the  closure  problem. 
The  use  of  transient  growth  methodology  in  transitional  flows  is  certainly  well-motivated.  However,  despite 
some  informative  results,  particularly  with  respect  to  the  match  between  the  most  amplified  spanwise  scales 
(occurring  for  streamwise-constant  disturbances)  and  experimental  observations  [del  Alamo  and  Jimenez, 
2006,  Hwang  and  Cossu,  2010],  the  analysis  of  infinitesimal  perturbations  growing  in  finite  time  in  a  system 
where  nonlinear  feedback  is  so  prevalent  struck  us  as  problematic  and  motivated  the  desire  to  look  beyond 
the  time-domain  behaviors  of  unforced  perturbations.  Turbulent  flows  experience  very  large  instantaneous 
deviations  from  any  known  fixed  point  solution  and  the  largest  growing  initial  condition  under  some  time 
horizon  will  surely  itself  be  perturbed  before  it  reaches  maturity.  We  suggest  that  things  look  simpler  in  the 
frequency  domain. 

Important  steps  were  made  in  this  direction  by,  among  others,  Bamieh  and  Dahleh  [2001]  and  Jovanovic  and 
co-authors  [Jovanovic  and  Bamieh,  2003,  2005].  There  is  a  link,  via  the  pseudospectrum,  between  the  non¬ 
normality  investigated  in  the  transient  growth  body  of  work  and  the  resolvent  norm,  which  we  will  describe 
further  in  what  follows  [Schmid  and  Henningson,  2001,  Trefethen  et  al.,  1993,  McKeon  and  Sharrna,  2010]. 
Those  studies  of  the  non-normality  of  the  linearized  Navier-Stokes  operator,  however,  did  not  go  so  far  as  to 
develop  a  full,  self-consistent  approach  to  wall-bounded  turbulence,  in  the  sense  that  links  to  experimental 
observations  were  not  made  and  the  connection  with  the  omitted  nonlinear  terms  was  not  established.  The 
resolvent  analysis  that  formed  a  central  part  of  the  work  performed  under  this  grant  seeks  to  fill  this  gap. 

Linear  models  of  pipe  flow,  e.g.,  Schmid  and  Henningson  [2001]  have  been  found  to  capture  general  char¬ 
acteristics  of  the  coherent  structures  present  in  turbulent  flow  but  are  unable  to  reproduce  the  change  in 
mean  flow  associated  with  transition  to  turbulence.  Linear  studies  have  shown  that  linear  non-normality  in 
subcritical  shear  flows  is  required  to  sustain  turbulence  [Henningson  and  Reddy,  1994]  and  so  are  the  terms 
linear  in  turbulence  fluctuations  [Lim  and  Kim,  2004].  The  non-normality  of  the  linear  Navier-Stokes  (LNS) 
operator  results  in  large  amplification  of  disturbances  and  transient  growth  of  initial  perturbations.  However, 
the  growth  of  the  most  amplified  structures  modify  the  mean  flow  in  a  way  that  reduces  the  amplification 
potential,  i.e. ,  the  non- normality.  Hence,  realistic  models  for  pipe  flow  transition  must  include  mechanisms 
leading  to  a  change  in  mean  flow  as  the  perturbations  develop. 

Since  the  blunting  of  the  velocity  profile  is  accurately  captured  by  the  full  Navier-Stokes  equations  (NSE), 
a  convenient  approach  is  to  go  back  to  the  NSE  and  invoke  simplifying  assumptions  in  order  to  reduce  the 
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complexity  of  the  model  instead  of  starting  from  the  LNS  and  adding  some  sort  of  nonlinearity.  Streamwise- 
constant  models  describe  the  evolution  of  the  three  components  of  velocity  in  a  plane  perpendicular  to 
the  mean  flow  and  are  equivalently  referred  to  as  2D/3C.  A  streamwise-constant  model  for  fully  developed 
(pipe)  flow  was  derived  by  Joseph  and  Tao  [1963]  and  shown  to  be  globally  stable  for  all  Reynolds  numbers. 
Thus,  the  2D/3C  model  has  the  useful  property  of  having  a  unique  fixed  point  corresponding  to  the  laminar 
flow.  A  stochastically  forced  2D/3C  model  formulated  in  terms  of  a  cross-stream  stream-  function  and  the 
deviation  of  the  streamwise  velocity  from  the  (linear)  laminar  profile  was  used  by  Gayme  et  al.  [2010]  to  study 
Couette  flow  and  successfully  captured  both  the  blunting  of  the  velocity  profile  and  structures  similar  to 
the  streamwise-elongated  vortices  and  streaks  observed  in  experiments.  In  general  terms,  the  stochastically 
forced  2D/3C  model  exploits  the  large  amplification  of  background  disturbances  due  to  the  non-normality  of 
the  linearized  operator  described  by  Farrell  and  Ioannou  [1998]  which  has  been  shown  to  reach  a  maximum  for 
streamwise-constant  disturbances  [Barnieh  and  Dahleh,  2001].  The  latter  authors  showed  that  streamwise- 
constant  disturbances  are  amplified  proportionally  to  Re 3  vs.  Re 3/2  for  streamwise-variant  disturbances. 
In  addition,  Jovanovic  and  Barnieh  [2005]  demonstrated  that  the  largest  amplification  is  only  obtained  by 
forcing  in  the  plane  perpendicular  to  the  mean  flow  and  is  observed  in  the  streamwise  velocity  component. 

Pipe  flow  is  well  suited  for  an  assumption  of  streamwise  invariance  since  streamwise-elongated  coherent  struc¬ 
tures  have  been  shown  to  play  an  important  role  during  transition  as  well  as  in  fully  developed  turbulence, 
e.g.,  Eckhardt  et  al.  [2007],  Kim  and  Adrian  [1999],  Morrison  et  al.  [2004].  The  streamwise-elongated  coher¬ 
ent  structures  in  pipe  flow,  both  in  the  near-wall  region  and  further  from  the  wall,  take  a  form  dominated  by 
quasi-streamwise  vortices  and  streaks  of  streamwise  velocity.  A  body  of  recent  work  in  the  literature  suggests 
a  connection  between  these  features  and  studies  of  the  LNS.  For  example,  the  most  (temporally)  amplified 
mode  of  the  LNS,  based  on  an  energy  norm,  is  streamwise-constant  with  an  azimuthal  wavenumber  n  =  1 
and  features  a  pair  of  counter-rotating  vortices  which  create  streaks  by  convecting  streamwise  momentum 
[Schmid  and  Henningson,  2001]. 

Traditionally,  the  later  stages  of  transition  to  turbulence  in  pipe  flow  have  been  characterized  by  the  creation 
of  puffs  and  slugs  [Wygnanski  and  Champagne,  1973].  Puffs  have  been  identified  as  the  flow  response  to 
large  amplitude  disturbances  at  low  Reynolds  number,  e.g.,  Re  ~  2000,  and  are  characterized  by  a  sharp 
trailing  edge  and  a  smooth  leading  edge  whereas  slugs  are  created  by  low  amplitude  disturbances  at  larger 
Reynolds  number,  Re  >  3000  and  have  sharp  leading  and  trailing  edges. 

3.2  Structure  of  the  turbulent  boundary  layer 

The  current  understanding  of  the  structural  and  statistical  nature  of  the  turbulent  boundary  layer  over  a 
large  range  of  Reynolds  number  rests  primarily  upon  spatial  information  from  direct  numerical  simulations 
(DNS)  and  particle  image  velocimetry  (PIV)  experiments,  and  temporal  information  from  experimental 
measurements  at  a  single  point,  yet  there  is  a  lack  of  simultaneous  temporal  and  spatial  measurements  to 
describe  the  time  evolution  of  the  structures  and  statistics  of  the  flow.  While  temporal  information  can  be 
extracted  from  DNS  data,  not  only  is  the  range  of  Reynolds  numbers  limited,  but  the  computational  resources 
required  to  analyze  a  time  resolved  flow  over  a  sufficiently  long  time  period  would  be  enormous.  Similarly, 
in  experiments,  most  commonly  used  instruments,  such  as  hot  wire  anemometers,  only  provide  temporal 
information  at  a  single  point,  where  the  recovery  of  spatial  information  requires  the  arduous  task  of  taking 
measurements  at  two  or  more  points  simultaneously  over  a  range  of  separations.  Planar  PIV  measurements 
provide  a  2D  spatial  representation  of  the  flow,  and  with  a  sufficiently  high  frame  rate,  as  shall  be  discussed, 
provide  simultaneous  spatial  and  temporal  measurements  allowing  a  statistical  and  structural  analysis  of  the 
turbulent  boundary  layer  over  two  spatial  dimensions  and  time. 

Often,  when  either  the  spatial  or  temporal  information  is  not  acquired,  a  conversion  between  the  two  is 
performed  using  Taylor’s  frozen  turbulence  hypothesis  [Taylor,  1938]  which  states  that  the  spatial  field  can 
be  reconstructed  from  the  temporal  field  if  the  convection  velocities  of  the  individual  eddies  or  scales  which 
compose  the  flow  are  known.  This  conversion  is  performed  assuming  the  eddies  are  “frozen”,  or  in  other 
words,  that  their  shape  does  not  evolve  significantly  over  the  distance  projected.  While  the  convection 
velocity  of  all  such  eddies  is  usually  assumed  to  be  equal  to  the  local  mean  velocity,  many  studies  have 
shown  that  this  is  not  always  a  valid  assumption  [Morrison  and  Kronauer,  1969,  Kim  and  Hussain,  1993, 
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Krogstad  et  al.,  1998,  Chung  and  McKeon,  2010,  del  Alamo  and  Jimenez,  2009,  LeHew  et  al.,  2010].  A 
recent  example  of  the  importance  of  knowing  the  correct  conversion  was  demonstrated  by  Monty  and  Chong 
[2009]  who  compare  channel  flow  computations  and  experiments  under  the  same  conditions  using  the  spatial 
spectrum  from  computations  and  the  converted  temporal  spectrum  from  experiments.  The  discrepancies 
between  the  spatial  and  temporal  spectra  were  put  into  better  agreement  using  a  scale  dependent  convection 
velocity  near  the  wall.  In  the  PIV  measurements  performed  under  this  grant,  information  was  recorded  in 
both  time  and  space,  so  no  conversion  was  necessary,  and  it  was  possible  to  test  the  validity  of  Taylor’s 
hypothesis  at  all  measurement  locations. 

The  use  of  time-resolved  PIV  to  address  the  need  for  both  spatial  and  temporal  measurements  in  the  turbulent 
boundary  layer  was  first  considered  by  Dennis  and  Nickels  [2008]  at  a  wall  normal  location  of  y/5  =  0.16, 
where  S  is  the  boundary  layer  thickness.  It  was  concluded  that  at  this  particular  wall  normal  location  with  a 
65  streamwise  held  of  view,  Taylor’s  hypothesis  holds,  although  deviation  was  noted  between  the  actual  and 
projected  velocity  fields  at  distances  beyond  4 5  downstream.  Based  on  the  works  cited  previously,  a  more 
notable  deviation  would  likely  appear  with  measurements  both  very  near  and  far  from  the  wall,  which  will 
be  addressed  in  the  current  study  over  a  larger  spatial  held  and  considering  the  entire  ( kx,kz,u> )  domain. 

Although  the  use  of  PIV  for  simultaneous  spatial  and  temporal  measurements  of  a  turbulent  boundary  layer 
is  promising,  there  are  a  number  of  issues  to  be  addressed.  First,  the  resolution  of  the  spatial  fluctuations  of 
the  how,  and  thus  the  resolution  of  the  spatial  spectrum,  is  limited  by  the  interrogation  window  size  which 
attenuates  small  scale  fluctuations  on  the  order  of  the  window  size  as  noted  by  Willert  and  Gharib  [1991]. 
A  technique  to  determine  the  optimal  window  size,  considering  both  the  attenuation  from  the  interrogation 
window  and  the  noise  introduced  by  the  PIV  algorithm,  was  proposed  by  Foucaut  et  al.  [2004].  An  outline 
of  all  of  the  effects  of  PIV  on  the  spatial  spectrum  including  methods  for  avoiding  spectral  leakage  and 
aliasing  in  space  is  presented  by  Tomkins  and  Adrian  [2005] .  The  study  of  the  space-time  correlation  in  hot 
and  cold  jet  hows  by  Wernet  [2007]  discusses  the  need  for  oversampling  and  low  pass  filtering  in  time  to 
avoid  temporal  aliasing  of  PIV  data.  Finally,  not  only  should  the  smallest  energetic  scales  in  the  how  be 
recovered,  but  also  the  largest  scales  found  to  extend  up  to  14  times  the  pipe  radius  in  experiments  by  Kim 
and  Adrian  [1999].  For  boundary  layer  measurements,  both  the  large  scale  motions  (LSMs)  with  an  energetic 
peak  at  Aa,  =  2  -  3 <5  and  the  superstructures,  with  an  energetic  peak  at  Xx  =  6<5  as  discussed  by  Hutchins 
and  Marusic  [2007],  Monty  et  al.  [2009],  and  Guala  et  al.  [2010],  where  Xx  is  the  streamwise  wavelength, 
must  be  resolved.  As  discussed  by  Balakumar  and  Adrian  [2007],  streamwise  scales  longer  than  3d  at  all 
wall  normal  locations  and  Reynolds  numbers  considered  contain  at  least  45%  of  the  turbulent  kinetic  energy 
in  the  zero  pressure  gradient  boundary  layer.  In  addition,  Ganapathisubramani  et  al.  [2003]  found  that 
hairpin  packets  of  comparable  streamwise  extent  contain  a  signihcant  amount  of  Reynolds  stress,  - uv ,  and 
thus  these  large  scales  must  be  resolved  to  accurately  represent  the  how.  As  noted  by  Monty  et  al.  [2009], 
while  the  6d  superstructure  peak  becomes  almost  undiscernable  outside  the  log  layer,  near  the  wall  this  peak 
is  prominent  and  structures  of  this  size  and  even  larger  must  be  resolved  to  fully  encompass  the  near  wall 
energy  and  dynamics. 

The  visualization  and  characterization  of  coherent  structures,  “organized  motions  that  are  persistent  in  time 
and  space  and  contribute  significantly  to  the  transport  of  heat,  mass,  and  momentum”  [Marusic  et  al.,  2010], 
in  turbulent  hows  is  a  powerful  tool  for  unwrapping  the  complicated  dynamics  involved  in  the  production 
and  sustenance  of  turbulence.  The  streaky  structures  observed  in  the  hydrogen-bubble  visualizations  of 
Kline  et  al.  [1967b]  in  the  turbulent  boundary  layer,  later  captured  in  simulations  by  Kim  et  al.  [1987],  led  to 
the  identification  of  the  autonomous  near-wall  cycle  analyzed  by  Jimenez  and  Pinelli  [1999]  as  well  as  Jeong 
and  Hussain  [1995]  and  Schoppa  and  Hussain  [2002],  recognized  as  the  major  producer  of  turbulent  kinetic 
energy  near  the  wall.  The  dynamics  of  turbulence  further  from  the  wall  are  often  associated  with  the  hairpin 
vortex  proposed  by  Theodorsen  [1952]  and  studied  in  the  smoke  visualizations  of  Head  and  Bandyopadhyay 
[1981].  A  comprehensive  review  of  the  hairpin  packet  paradigm  is  given  by  Adrian  [2007].  The  literature 
describing  coherent  structure  is  broad,  and  we  review  only  the  papers  most  relevant  to  the  current  study 
here. 

The  advent  of  direct  numerical  simulation  (DNS)  and  experimental  techniques  such  as  particle  image  ve- 
locimetry  (PIV)  have  allowed  not  only  a  qualitative,  but  also  a  quantitative  measure  of  the  flow  field  as¬ 
sociated  with  these  coherent  structures.  Adrian  et  al.  [2000]  found  that  hairpin  vortices  tend  to  form  into 
streamwise  aligned  packets  using  PIV  in  wall-normal  planes,  and  later,  Tomkins  and  Adrian  [2003]  analyzed 
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these  structures  in  wall-parallel  planes.  Work  by  Ganapathisubramani  et  al.  [2003]  showed  that  not  only 
were  hairpin  packet  structures  prevalent  in  the  flow,  but  they  also  produced  significant  amounts  of  Reynolds 
stress,  and  were  thus  linked  with  the  production  of  turbulence. 

In  depth  analysis  of  the  size  and  organization  of  these  hairpin  vortices  was  performed  by  Carlier  and  Stanislas 
[2005]  using  planar  2D2C  (two-dimensional  velocity  field  with  two  in-plane  components  of  velocity)  PIV  in 
wall-normal,  wall-parallel,  and  upstream-  and  downstream-tilted  (135°  and  45°  to  the  free-stream,  respec¬ 
tively)  cross-stream  planes  and  stereo-PIV  in  a  vertical  cross-stream  plane  of  the  turbulent  boundary  layer 
over  a  range  Reg  =  7500-19000,  where  Reg  =  OU^Iv  is  the  momentum  thickness  Reynolds  number,  6  is  the 
momentum  thickness,  Uoo  is  the  free-stream  velocity,  and  v  is  the  kinematic  viscosity.  It  was  found  that 
vortices  have  a  radius  ranging  from  20-30  and  that  their  properties  scale  universally  in  inner  units  in  the 
log  region  for  moderate  Reynolds  number  flows,  where  8V  -  vluT  is  the  viscous  length  scale  and  uT  -  \J tw/ p 
is  the  friction  velocity  where  tw  is  the  shear  stress  at  the  wall  and  p  is  the  fluid  density.  Additionally,  stereo- 
PIV  measurements  from  the  same  authors  indicated  that  vortex  structures  identified  in  the  cross-stream 
plane  travel  near  the  local  mean  velocity  on  average. 

Hutchins  et  al.  [2005]  also  performed  stereo-PIV  measurements  in  upstream  and  downstream  tilted  planes, 
135°  and  45°  to  the  free-stream,  respectively.  In  the  upstream  tilted  plane,  trains  of  counter  rotating  vortex 
structures  were  found  where  the  orientation  of  pairs  producing  an  ejection  was  preferred,  consistent  with 
the  hairpin  packets  proposed  by  Adrian  et  al.  [2000].  Furthermore,  the  spacing  of  vortex  pairs  scaled  in 
outer  units  and,  using  the  streamwise  velocity  correlation  coefficient,  structures  were  inferred  to  uncouple 
or  detach  from  the  wall  once  they  grew  beyond  the  log-layer,  and  thus  the  detachment  point  scaled  in 
outer  units  as  well.  Finally,  the  dual-plane  stereo  PIV  experiments  of  Ganapathisubramani  et  al.  [2006] 
helped  to  elucidate  the  orientation  of  vortex  structures  by  obtaining  the  complete  vorticity  vector.  These 
measurements  verified  that  the  most  frequently  observed  vortex  inclination  angle  ranged  from  42°  to  46°, 
showing  the  predominance  of  forward  or  downstream-leaning  vortices  near  the  wall  with  a  A  shaped  vortex 
being  the  preferred  representation,  and  demonstrated  that  the  density  of  backward  leaning  cores  is  nearly 
constant  with  y,  while  the  number  of  forward  leaning  cores  decreases  with  y. 

While  much  progress  has  been  made  by  obtaining  static  spatial  footprints  of  coherent  structures  in  the 
turbulent  boundary  layer,  there  has  been  little  study  of  the  evolution  of  these  coherent  structures  in  time. 
This  is  difficult  experimentally  due  to  the  need  for  a  full-field  velocity  measurement  in  a  large  field  of 
view,  along  with  fine  temporal  resolution;  usually  only  one  of  the  two  is  possible  with  a  given  experimental 
technique.  While  temporal  resolution  may  not  be  an  issue  in  DNS,  the  memory  requirements  to  store  a 
significantly  long  time  history  to  track  coherent  structures  and  recover  a  statistically  significant  measure 
would  be  large. 

Many  researchers  have  identified  vortical  structures  or  their  signature  in  the  Reynolds  stress  in  direct  numer¬ 
ical  simulations  and  explored  the  relationships  between  them,  e.g.  Jimenez  and  co-workers,  del  Alamo  et  al. 
[2006],  Lozano-Duran  et  al.  [2012],  albeit  constrained  to  spatial  averaging  rather  than  temporal  tracking. 
Here  we  focus  on  temporally-resolved  measurements,  which  are  predominantly  experimental  at  present. 

One  alternative  to  the  stringent  resolution  and  field  of  view  requirements  in  experiment  is  to  acquire  two-  or 
three-velocity-component,  two-dimensional  fields  and  extrapolate  into  three  dimensions.  In  order  to  better 
understand  the  3D  structure  and  evolution  of  vortex  structures,  Dennis  and  Nickels  [2011b]  and  Dennis  and 
Nickels  [2011a]  used  stereo  PIV  in  the  cross-stream  plane  of  a  turbulent  boundary  layer  to  obtain  “pseudo-3D” 
velocity  fields  and  reconstruct  the  spatial  development  of  hairpin  vortices  and  large-scale  coherent  structure 
using  Taylor’s  hypothesis.  This  hypothesis  is  known  to  be  accurate  over  a  large  portion  of  the  boundary 
layer,  thus  providing  a  fairly  accurate  representation  of  the  3D  structure.  By  conditioning  the  velocity  field 
on  a  swirl  threshold  at  a  particular  height,  the  ensemble  relationship  between  structures  was  developed. 

Recent  time-resolved  tomographic  PIV  measurements  in  a  1.8  x  1.8  x  0.19  5  (streamwise,  spanwise,  wall- 
normal)  field  of  view  by  Elsinga  and  Marusic  [2010]  have  satisfied  both  the  spatial  and  temporal  resolution 
requirement  and  allowed  the  measurement  of  the  average  lifetime  of  coherent  structures,  which  they  measure 
to  be  14.3(5/C/oo  based  on  the  period  of  oscillation  between  the  four  basic  flow  topologies  (inferred  from 
the  invariants  of  the  local  velocity  gradient  tensor;  see  Chong  et  al.  [1990]),  where  S  is  the  boundary  layer 
thickness.  Using  the  same  data  set,  Elsinga  et  al.  [2012]  were  also  able  to  track  vortex  structures  in  time, 
and  observed  non-negligible  wall-normal  displacement  of  the  structures  during  a  typical  trajectory. 
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4  Critical  layer  framework 


In  this  section,  a  review  of  the  critical  layer/resolvent  framework  results  obtained  under  this  grant  is  given. 
This  encompasses  material  published  in  three  papers  in  the  Journal  of  Fluid  Mechanics  [McKeon  and  Sharma, 
2010,  Sharma  and  McKeon,  2013b,  Moarref  et  al.,  2013].  The  initial  discussion  closely  follows  the  set  up 
in  the  (invited)  Physics  of  Fluids  review  paper  McKeon  et  al.  [2013].  The  origins  of  this  analysis  partly  lie 
in  nonlinear  robust  control  theory,  but  a  differentiating  feature  is  the  connection  with,  and  prediction  of, 
state-of-the-art  understanding  of  velocity  statistics  and  coherent  structures  observed  in  real,  high  Reynolds 
number  flows. 


4.1  Turbulence  as  a  directional  amplifier 

The  full  analysis  described  in  McKeon  and  Sharma  [2010]  was  performed  for  pipe  flow  primarily  because 
of  the  geometrical  simplicity,  the  availability  of  high  Reynolds  number  experimental  data  for  comparison  of 
statistical  results  [McKeon  et  al.,  2004,  Morrison  et  al.,  2004]  and  the  potential  for  extrapolation  to  other 
canonical  flows.  Turbulent  flow  through  pipes  is,  of  course,  important  for  the  transport  of  fluids  such  as 
oil  and  natural  gas,  numerous  natural  and  biomedical  applications,  and  is  also  highly  relevant  to  the  study 
of  other  canonical  flows.  Transition  to  turbulence  in  a  pipe  is  still  not  completely  understood  Hof  et  al. 
[2010],  but  the  pipe  offers  the  analytical  benefits  of  statistical  homogeneity  in  the  streamwise  direction  and 
a  simple  constraint  on  the  azimuthal  wavenumber.  The  coordinate  system  in  a  long,  straight  pipe  of  circular 
cross-subsection  for  the  following  analysis  is  shown  in  figure  12. 

The  conceptual  picture  underlying  our  analysis  is  that  of  a  nonlinearity  “feeding  back”  to  excite  the  linear 
wave  propagation  dynamics,  which  in  turn  drive  the  nonlinearity.  As  such,  we  concentrate  on  the  linear 
amplification  aspect  of  the  turbulent  process  and  examine  the  corresponding  linear  transfer  function,  con¬ 
sidering  its  action  at  each  wavenumber/frequency  combination  (described  below),  which  is  shown  to  provide 
highly-directional  amplification.  This  picture  is  closely  related  to,  but  distinct  from,  the  sector-bounding 
analysis  of  the  NSE  for  control  [Sharma  et  al.,  2011]  and  model  reduction  [Sharma,  2009]. 

A  full  solution  of  the  NSE  would  predict  the  turbulent  mean  velocity  profile.  However,  a  significant  restriction 
on  the  scope  of  our  modeling  efforts  is  obtained  by  assuming  the  mean  profile  is  known  and  concentrating  on 
the  action  of  the  resolvent  operators  in  which  it  appears.  With  the  particular  choice  of  pipe  flow,  experimental 
measurements  from  the  Princeton  Superpipe  are  available  up  to  conditions  representative  of  high  Reynolds 
number  applications  [McKeon  et  al.,  2004]  and  the  use  of  a  turbulent  eddy  viscosity  can  be  avoided. 

Pipe  flow  affords  the  simplification  of  a  priori  knowledge  that  the  optimal  bases  for  modeling  the  flow  field 
in  the  homogeneous  directions  are  the  eigenmodes,  which  are  the  Fourier  modes,  such  that  the  flow  can  be 
represented  as  a  summation  of  finite-amplitude  traveling  wave  fluctuations  over  a  range  of  two-dimensional 
spatial  wavenumbers  and  temporal  frequencies  about  the  turbulent  mean  profile,  with  a  natural  connection 
to  the  gain  analysis  described  previously.  Another  obvious  attraction  of  Fourier  decomposition  is  the  direct 
connection  to  the  spectral  measurements  typically  made  in  real  flows.  At  its  heart,  the  analysis  that  follows 
seeks  to  establish  the  optimal  basis  for  decomposition  in  the  inhomogeneous  wall-normal  direction,  or  the 
radial  coherence  associated  with  particular  traveling  waves.  The  analysis  has  been  compared  to  proper 
orthogonal  decomposition,  which  is  a  well-known  technique  for  decomposing  velocity  field  data  to  optimally 
capture  the  energy.  While  there  is  a  superficial  correspondence,  our  analysis  differs  in  that  the  decomposition 
is  predictive,  since  it  is  applied  to  the  NSE  and  not  pre-existing  data. 


4.1.1  Formulation  for  the  analysis 

Our  analysis  begins  with  the  non-dimensional  NSE  for  fully-developed,  incompressible  pipe, 

<9*u  =  -S/p  -  u  •  Vu  +  Re-1  V2u,  V-u=0  (1) 

We  retain  the  boundary  layer  terminology  by  fixing  y  =  1  -  r,  and  u ,  v(=  —v')  and  w  as  corresponding  to 
the  streamwise,  wall-normal  and  azimuthal  velocities  such  that  u  =  ( v',w,u ),  as  shown  in  figure  12.  The 
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Figure  12:  A  schematic  of  pipe  geometry  and  nomenclature. 


Reynolds  number  is  defined  as 


n  f^bulk-D 

Re  = - , 

v 

where  v  is  a  constant  viscosity,  D  the  pipe  diameter  (equal  to  2 R)  and  the  bulk,  volume-averaged  velocity 
is  U bulk •  We  also  introduce  the  Karman  number, 


R+  = 


DuT 

~2v~' 


Here,  uT  =  yfr^Jp  is  the  friction  velocity,  tw  is  the  mean  wall  shear  stress  and  p  is  the  density. 

The  NSE  are  invariant  under  translation  in  t  (time),  x  and  6.  In  view  of  these  symmetries,  we  Fourier 
transform  in  time  and  the  two  homogeneous  spatial  directions,  introducing  for  convenience  the  wavenum- 
ber/wavespeed  triplet  I\  =  ( k,n,c )  with  k  =  k'R,  n  =  n'R  and  wavespeed  c  =  ui/k,  with  ui  =  ui'R/ uo\y+=n+ 
(where  the  prime  denotes  the  dimensional  variables  and  u0|y+=fl+  is  the  centerline  velocity),  and  let  K  •  x  = 
kx  +  nO  -  uit  (notationally  convenient).  In  addition,  we  specially  define  Kq  as  ( k,n,u> )  =  (0,0,0)  and  note 
that  Uk0  (u0)  is  the  turbulent  velocity  field  averaged  over  space  and  time.  Thus  K  with  all  positive  elements 
refers  to  a  downstream  traveling  wave  (helical  in  the  cylindrical  pipe  geometry).  As  yet,  only  the  wall-normal 
direction  remains  untransformed;  the  problem  is  to  find  a  suitable  basis  for  these  functions  of  r.  The  velocity 
field  (and  forcing  field)  is  expressed  as  a  sum  of  harmonic,  radially  varying  traveling  waves. 

Xoo  r*  oo 

/  uK(r)ezKxdkduj  (2) 

OO  J  —  oo 


We  define  u  and  f  as 


u  =  u  -  u0  (3) 

f  =  -u  ■  Vu  (4) 

with  u k(x)  and  fA'(?’)  the  (radially  varying)  Fourier  coefficients  thereof. 

The  NSE  assuming  fully  developed  flow  (i.e.  real  elements  of  K)  can  then  be  written  in  a  divergence-free 
basis  as 


-iuj\iK  =  Ckvlr  +  f/c,  (k,n,oj)  *  (0,0,0)  (5) 

0  =  f0  -  u0  •  Vu0  +  i?e_1V2uo  (6) 

Note  that  in  this  formulation,  the  linearized  Navier-Stokes  operator,  Ck ,  naturally  involves  the  turbulent 
mean;  yet  we  have  performed  no  (linearized)  stability  analysis  around  a  fixed  point.  Thus  the  u/<-s  are 
exactly  the  Fourier  coefficients  of  the  turbulent  fluctuations  relative  to  the  mean  profile  at  all  y.  Previous 
work  concerning  analysis  of  Ck  and  globally  optimal  perturbations  has  typically  employed  a  turbulent 
viscosity  to  account  for  the  turbulent  transport  that  would  be  excluded  by  a  linearization  procedure.  In 
this  analysis,  the  treatment  of  f  explicitly  accounts  for  this  phenomenon,  so  the  resolvent  contains  only  the 
turbulent  mean  velocity  profile  itself. 
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Figure  13:  A  block  diagram  showing  the  network  of  resolvents,  Hk ,  acting  on  radially- varying  traveling 
waves  of  different  wavenumber  and  frequency.  An  input-output  relationship  can  be  written  for  all  K,  but 
ux  +  0  for  a  finite  range  of  K  in  real  flows.  FT  and  IFT  denote  the  Fourier  transform  and  inverse  Fourier 
transform,  respectively,  uq  and  /o  describe  the  equation  for  the  turbulent  mean  profile. 

A  trivial  manipulation  of  equation  5  gives  the  response  of  the  flow  at  a  particular  wavenumber  combination 
subjected  to  harmonic  forcing  f/<  arising  from  the  interaction  between  other  wavenumbers, 

uk  =  Hk^k  =  (-iwl  -  Ck) (7) 

where  Hk  is  known  in  the  literature  as  the  resolvent  [Schmid,  2007].  Using  the  form  for  Ck  given  in  Meseguer 
and  Trefethen  [2003],  we  can  write 

i(kuo  -  gj)  -  Re~xD  -2inr~2  Re-1  0 

Hk  =  2inr~2Re~1  i(kuo  -  u>)  -  Re~1D  0  (8) 

-dru o  0  z(fcu0  -  w)  -  Re~1{D  +  r~2) 

with  D  =  d2  +  r~1dr-r~2(  n2  + 1)  -&r,  and  the  states  being  the  radial,  azimuthal  and  axial  velocities  expressed 
in  a  divergence-free  basis.  Thus  equations  5  and  6  can  be  visualized  as  the  block  diagram  in  figure  13,  forming 
a  complete  representation  of  the  NSE  in  this  basis.  The  resolvents,  Hk,  can  be  derived  from  the  NSE  and 
are  interconnected  through  the  nonlinear  terms.  Note  that  (as  stated  in  the  caption)  while  an  input-output 
relationship  can  be  written  for  all  K,  U/^-  +  0  for  a  finite  range  of  K  in  real  flows. 

Although  the  resolvents  are  linear ,  the  nonlinear  forcing  internal  to  the  system  is  explicitly  retained  so  this 
is  not  a  linearized  analysis.  Additional  effects  of  the  nonlinearity  are  bound  up  in  the  use  of  the  turbulent 
mean  velocity  profile  in  the  resolvent.  Thus  the  analysis  is  fundamentally  different  from  linear  stability 
analysis:  it  concerns  the  response  of  a  forced  system  (in  reality  lightly  forced  because  the  leading  singular 
values  in  general  give  rise  to  large  amplification),  where  the  forcing  is  required  to  sustain  a  response  in  a 
linearly  stable  system.  Thus,  concepts  relevant  to  the  study  of  disturbances  in  inviscid,  linearized  laminar 
flow  can  be  extended  to  the  turbulent  case,  with  the  understanding  that  in  the  latter  case  the  waves  are 
lightly  damped  and  would  asymptotically  decay  in  the  absence  of  forcing  f. 

Figure  13  describes  an  interconnected  series  of  linear  sub-systems  with  resolvents  of  known  form.  Thus  the 
formulation  of  equations  5  and  6  lends  itself  naturally  to  separate  analysis  of  the  individual  sub-systems 
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Figure  14:  An  illustration  of  the  effect  of  matrix  M  on  the  unit  circle  (left  plot  mapped  to  centre  plot). 
The  singular  values  are  the  radii  of  the  resulting  ellipse  and  the  singular  vectors  give  the  rotation  of  the 
ellipse.  a;  are  the  columns  of  A  and  b;  are  the  columns  of  B.  The  rightmost  plot  shows  the  effect  of  the 
optimal  rank-1  approximation  M  to  the  mapping  M  on  the  unit  circle. 


associated  with  each  K,  under  the  constraint  that  the  correct  amplitudes  for  each  ux  are  required  to 
sustain  the  (assumed)  mean  velocity  profile  in  the  fully-connected  system.  Without  satisfying  this  constraint 
associated  with  the  excitation  of  all  u^-  through  the  nonlinear  term,  progress  can  still  be  made  in  terms  of 
determining  the  radial  distributions  of  velocity  for  each  K,  since  the  linearity  of  the  system  allows  for  simple 
amplitude  scaling  of  the  response  mode  shapes.  In  effect,  then,  we  can  consider  unstructured  (unit)  forcing 
of  each  sub-system,  in  knowledge  that  the  output  response  modes  can  be  linearly  superposed  to  return  to 
a  representation  of  the  full  velocity  field  subject  to  the  constraint  that  the  response  mode  amplitudes  must 
ultimately  satisfy  figure  13.  This  is  another  significant  advantage  of  our  approach:  the  ability  to  dissect 
the  flow  into  manageable  “lumps”  (linear  sub-systems)  that  can  be  analyzed  in  isolation  and  then  simply 
reassembled.  Further,  this  formulation  proves  amenable  to  approximation  of  the  resolvent  at  each  K  that 
leads  to  a  significant  reduction  in  complexity  using  standard  linear  systems  techniques. 


4.1.2  Optimal  low-rank  approximation  of  the  resolvent  for  inhomogeneous  coordinates 

The  symmetries  associated  with  the  pipe  geometry  imply  that  the  resolvent  is  normal  under  the  integrals 
over  the  wall-parallel  coordinates  and  time  (see  Sharma  and  McKeon  [2013a]  for  further  discussion  on  this 
point).  In  these  homogeneous  directions,  then,  the  forcing  and  response  modes  are  equal,  orthogonal  and 
are  the  Fourier  modes,  which  then  represent  the  appropriate  choice  for  basis  functions.  However,  this  is  not 
true  in  the  wall-normal  direction:  the  presence  of  the  wall  leads  to  the  loss  of  symmetry  with  the  implication 
that  the  forcing  and  response  modes  are  no  longer  equal  and  the  fluctuations  may  now  gain  energy  from  the 
interaction  with  the  mean  flow.  The  potential  for  momentum  production  due  to  this  interaction  is  quantified 
by  this  loss  of  orthogonality.  Since  the  resolvent  is  not  normal  with  respect  to  the  integral  in  this  direction, 
the  Fourier  basis  is  no  longer  the  optimal  choice  under  the  energy  norm.  We  seek  a  basis  that  is  optimal  with 
respect  to  the  magnitude  of  response  to  forcing  in  order  to  continue  our  gain-based  analysis.  To  find  this 
basis,  we  use  the  SVD.  For  the  reader  less  familiar  with  the  SVD,  we  give  a  brief  example  of  its  underlying 
action  and  interpretation  in  the  next  two  paragraphs  which  the  more  accustomed  reader  may  safely  skip. 

The  singular  value  decomposition  (SVD)  is  a  well-known  matrix  decomposition  that  splits  any  matrix  M 
into  the  product  M  =  ASB*  where  A  is  a  unitary  matrix  (AA*  =  A*A  =  I),  S  is  a  diagonal  real  matrix, 
and  B  is  the  conjugate  transpose  of  another  unitary  matrix,  all  of  appropriate  dimension,  as  in  the  example 
of  figure  14.  The  diagonal  elements  of  S,  cq,  are  ordered  such  that  a,  >  cq+i  and  are  called  the  singular 
values.  When  applied  to  statistical  data,  the  SVD  is  commonly  known  as  principal  component  analysis, 
Karhunen-Loeve  decomposition,  or  proper  orthogonal  decomposition  (POD).  Comparison  across  application 
domains  is  not  very  meaningful  however,  since  the  object  of  optimization  is  quite  different. 

The  singular  values  contained  in  the  diagonal  matrix  S  quantify  the  underlying  dimensionality  of  the  mapping 
being  decomposed.  Where  cr2  is  small  compared  to  a i,  as  in  the  simple  example  in  figure  14,  the  full  mapping 
of  a  unit  circle  at  input  to  an  ellipse  at  output  is  well  approximated  by  the  lower-dimensional  mapping  to 
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a  line.  This  example  illustrates  the  use  of  the  SVD  to  find  lower-rank  approximations  to  linear  mappings 
between  higher-dimensional  spaces.  The  SVD  generalizes  to  apply  to  linear  operators  on  Hilbert  spaces 
[Young,  1988]  where  it  becomes  known  as  the  Schmidt  decomposition. 

Hopefully  by  now  the  application  of  the  SVD  in  our  context  becomes  obvious:  since  we  hold  that  turbulence 
can  be  captured  by  a  low-order  amplifying  process,  we  must  find  the  “best”  low-rank  approximation  to  the 
resolvent,  where  what  is  meant  by  “best”  is  quantifiable  in  a  defined  sense.  The  result  is  simply  a  truncation 
of  the  Schmidt  decomposition  of  the  resolvent,  or  the  identification  of  the  right  singular  vectors  (inputs) 
that  give  rise  to  the  most  amplified  left  singular  vectors  (outputs).  The  relative  amplification  factor  for  unit 
input  amplitude  is  given  by  the  associated  singular  value.  We  term  the  left  and  right  singular  vectors  the 
velocity  response  and  forcing  modes,  respectively;  the  first  singular  forcing  modes  are  the  most  significant 
in  that  they  can  be  expected  to  dominate  the  full  observed  response  at  a  given  K  if  the  resolvent  is  truly 
low  rank.  We  confirm  this  assumption  in  the  results  shown  below,  noting  that  the  true  output  from  each 
resolvent  will  be  the  product  of  forcing  amplitude  and  singular  value. 

The  SVD  of  the  resolvent  before  Fourier  decomposition  is  imposed  would  itself  naturally  result  in  a  de¬ 
composition  into  the  (optimal)  Fourier  basis  in  the  symmetric  directions  for  the  reasons  discussed  earlier. 
For  simplicity,  we  take  that  stage  as  given  and  consider  the  optimal  low  rank  approximation  only  in  the 
inhomogeneous,  wall-normal  direction. 

The  SVD  of  the  resolvent  for  a  particular  I\  can  be  written  in  terms  of  the  forcing  (<p)  and  response  (ip) 
modes  associated  with  the  response  of  the  flow  as  follows: 

(-iwl  -  CkY1  =  YipiK(r)alK(pm(r)  (9) 

1  =  1 


with  the  orthogonality  condition 

( <PlK(r),<pmK(r))r  =  Slm ,  (lplK(r),1pmK(r))r  =  Sim,  (10) 

where  (v)r  indicates  the  inner  product  over  r,  such  that  the  sets  of  both  forcing  and  response  modes  are 
normalized  with  respect  to  the  energy  integrated  over  the  radius  of  the  pipe  and  are  orthogonal.  The  singular 
values  are  sorted,  with  aitK  >  cti+i,k  >  0.  For  stable  Ck ,  this  decomposition  exists  for  real  cj.  It  is  also  unique 
up  to  a  pre-multiplying  unitary  complex  factor  on  both  bases  corresponding  to  a  phase  shift  and  up  to  the 
ordering  in  l  of  aw s,  hence  we  fix  the  relative  phases  with  respect  to  the  first  coefficient. 

The  basis  pair  defined  by  (pm  and  ipw  can  be  used  to  decompose  arbitrary  forcing  and  the  resulting  velocity 
at  any  particular  frequency/wavenumber  component  K ,  such  that 

M?’)  =  Y,<PlK(r)xiK 
1  =  1 
oo 

uk(t)  =  Y,  criK^iK(r)xiKi 
i=i 

where  the  y’s  correspond  to  the  coefficients  of  projection  of  the  forcing  onto  the  forcing  modes. 

The  radial  shape  of  harmonic  forcing  that  gives  the  largest  disturbance  energy  (given  unit  forcing  amplitude) 
at  a  particular  K  is  f/<-  =  <Pik,  with  a  gain  of  am-  The  next  largest  arises  from  f#-  =  (p2K  and  so  on,  at  a 
particular  wavenumber  pair  and  frequency.  The  corresponding  flow  response  modes  are  given  by  the  related 
u/<-  =  ipiK,ip2K,  etc.  This  decomposition  of  general  input  forcing  is  illustrated  in  figures  15(a)  and  (b). 

At  this  stage,  the  analysis  is  complete  in  the  divergence-free  basis,  in  the  sense  that  no  modeling  assumptions 
have  been  made,  other  than  that  the  velocity  field  is  statistically  homogeneous  in  x ,  6  and  t.  Linearity  of 
the  system  implies  that  the  resolvent  can  be  analyzed  under  unstructured  forcing,  i.e.  isolating  individual 
resolvents  as  shown  in  figure  15(a),  but  the  appropriate  amplitudes  of  the  velocity  response  modes  must  be 
determined  before  the  connected,  self-sustaining  representation  of  figure  13  can  be  reassembled  (the  topic 
of  ongoing  work).  The  optimal  low-rank  approximation  to  the  resolvent  is  defined  to  arbitrary  accuracy 
related  to  the  number  of  singular  modes  included.  There  is  a  further  technical  assumption  that  the  resolvent 
is  compact.  The  simplest  approximation  is  rank-1,  in  which  the  resolvent  is  modeled  using  only  the  first 
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Figure  15:  Block  diagram  for  the  equations  (11)  and  (12)  referring  to  the  linear  equations  projected  onto 
the  forcing  and  response  modes.  Diagram  (a)  depicts  the  ‘true’  resolvent  and  (b)  depicts  the  Schmidt 
decomposition  of  the  resolvent.  The  rank-1  approximation  is  depicted  in  (c).  This  approximation  can  be 
performed  for  each  lower  sub-system  in  figure  13. 


singular  modes  (figure  15(c));  this  approximation  is  likely  to  be  good  if  cr,  K  »  alK  for  all?  >  1.  While  the 
assumption  that  a  rank-1  approximation  is  sufficient  to  model  wall  turbulence  is  quite  strong  to  begin  with, 
it  is  one  of  the  key  outcomes  of  this  formulation  that  it  yields  surprisingly  good  agreement  with  observations 
from  real  flows  for  the  vast  majority  of  observed  K  combinations.  Of  course,  a  higher-rank  approximation 
may  simply  be  made  by  the  same  process,  should  higher  fidelity  be  required. 

4.1.3  Response  of  radially- varying  traveling  waves  and  the  connection  with  critical  layer 
theory 

The  singular  response  modes  were  computed  using  a  MATLAB  implementation  of  the  decomposition  de¬ 
scribed  in  equations  in  the  previous  subsection  and  the  numerical  scheme  and  divergence-free  basis  for  pipe 
flow  proposed  by  Meseguer  and  Trefethen  [2003],  modified  for  a  turbulent  mean  flow  as  in  McKeon  and 
Sharma  [2010].  Full  details  of  the  computational  approach  can  be  found  in  these  two  papers  and,  for  brevity, 
are  not  repeated  here.  The  computational  expense  is  exceedingly  low,  even  for  analysis  of  Reynolds  number 
up  to  R+  =  1.97  x  104,  the  highest  Re  study  we  have  performed  to  date. 

McKeon  and  Sharma  [2010]  investigated  the  variation  of  the  characteristics  of  the  radially-varying  traveling 
waves  over  a  range  of  K  and  made  some  observations  with  respect  to  variation  with  singular  value.  Recall 
that  the  cr;  gives  the  amplification,  and  complex  functions  <j>i  and  ipi  give  the  wall-normal  variation  of  the 
amplitude  and  phase  of  the  forcing  and  response  mode,  respectively  (where  the  phase  of  the  response  mode 
is  relative  to  the  forcing).  In  particular,  the  orthogonality  requirements  given  in  equation  10  require  an 
increase  in  the  number  of  maxima  in  the  radial  amplitude  function  u k  (r)  with  increasing  l  (with  analogy  to 
the  orthogonality  of  Fourier  modes),  as  shown  in  the  comparison  of  velocity  response  modes  with  increasing 
l  in  figure  2  of  McKeon  and  Sharma  [2010].  (An  example  of  the  variation  of  singular  values  with  l  is  shown 
in  figure  3  of  that  paper.)  The  combination  of  the  increasing  complexity  of  such  structures  and  the  rapid 
drop-off  in  singular  value  observed  with  increasing  l  support  the  rank-1  hypothesis  for  a  wide  range  of  K ;  in 
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Figure  16:  (a)  Variation  of  the  location  of  the  peak  streamwise  velocity  with  increasing  wavespeed,  c,  for 
the  first  velocity  response  mode  with  ( k,n )  =  (1,10)  at  R+  =  1800  (solid  line).  The  turbulent  mean  profile 
normalized  by  the  centerline  velocity,  Uo/uo|y+=/j+ ,  at  this  Reynolds  number  is  also  shown  in  dotted  gray  for 
reference,  (b)  The  corresponding  first  singular  values. 


the  rest  of  this  paper,  we  consider  only  the  output  of  the  rank-1  model  of  the  NSE  at  each  K  value  (figure 
15(c)),  which  we  call  the  first  velocity  response  modes. 

We  summarize  the  general  features  of  the  first  velocity  response  modes  in  figures  16-17,  for  a  particular 
(fc,n)  =  (1, 10)  at  R+  -  1800.  Figure  16(a)  shows  the  wall-normal  variation  of  the  location  of  the  magnitude 
of  the  streamwise  velocity  perturbation  with  increasing  wavespeed  c,  compared  with  the  local  turbulent  mean 
velocity  profile  Uq.  The  variation  with  c  of  the  associated  first  singular  value,  er i  is  shown  in  figure  16(b); 
note  that  for  this  K  (and  a  wide  range  of  I\  in  general),  the  magnitude  of  the  first  singular  value  is  extremely 
large,  meaning  that  only  very  small  forcing  amplitude  will  be  required  to  observe  a  strong  velocity  response. 

Three  different  types  of  modes  can  be  identified  in  figure  16(a).  The  velocity  response  mode  shapes  u k 
corresponding  to  these  three  classes  of  velocity  response  modes,  identified  here  by  (i)-(iii),  are  shown  in 
figure  17.  The  variation  of  the  azimuthal  component  is  not  shown  since  it  is  simply  determined  from  the 
summation  of  left-  and  right-going  propagating  waves  described  earlier  and  the  continuity  equation,  and  is 
out  of  phase  with  u  in  both  streamwise  and  azimuthal  directions.  When  c  is  small  in  region  (i),  the  peak 
magnitude  in  all  components  stays  at  an  approximately  constant  wall- normal  distance,  y+  =  yuTlv.  We 
term  these  “attached  modes”,  because  the  footprints  of  the  modes  reach  down  to  the  wall.  For  large  c  in 
region  (iii),  the  location  of  the  peak  amplitude  tracks  the  local  mean  velocity,  i.e.  c  =  Uo  (y)  and  there  exists 
a  critical  layer.  The  mode  is  localized  around  the  critical  layer,  such  that  these  “critical  modes”  rapidly 
become  detached  from  the  wall  with  increasing  c.  Some  special  significance  will  be  given  to  the  class  of 
attached,  critical  modes  (ii),  i.e.  the  slowest  converting  modes  that  achieve  critical  status. 

Note  also  that  the  Reynolds  stress  associated  with  the  velocity  response  modes  exhibits  slightly  different 
behavior,  lifting  from  the  wall  to  provide  a  radially-localized  distribution  with  contributions  at  Kq  and  2 K. 

We  will  not  comment  here  on  the  first  singular  forcing  mode  shapes  corresponding  to  the  velocity  response 
modes  shown  in  figure  17,  beyond  saying  that  the  forcing  occurs  locally  to  the  velocity  response  in  a  radial 
sense,  and  that  understanding  which  pairs  of  other  Ks  supply  a  component  of  forcing  in  these  directions  is 
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Figure  17:  Representative  velocity  response  mode  shapes  for  the  “attached”,  “attached  and  critical”,  and 
“critical  modes”,  shown  for  conditions  corresponding  to  points  (i)  -  (Hi)  in  figure  16,  i.e.  (k,n)  =  (1,10), 
c-  0.25,0.67  and  0.85,  and  R+  =  1800.  Top:  u(y+ );  middle:  v(y+);  bottom:  uv(y+ ). 
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equivalent  to  understanding  how  to  close  the  block  diagram  of  figure  13. 

The  physics  underlying  these  distinct  classes  of  response  modes  can  be  identified  by  returning  to  the  full  rank 
resolvent  of  equations  7  and  8  (and  noting  that  the  scaling  phenomena  observed  are  reflected  in  the  rank-1 
approximation).  There  are  two  distinct  mechanisms  by  which  the  resolvent  can  lead  to  large  amplification 
[McKeon  and  Sharma,  2010]:  the  non-normality  associated  with  the  u9rUo  term  coupling  the  radial  pertur¬ 
bation  to  the  mean  shear  (requiring  a  mode  covering  a  wider  range  in  the  wall-normal  direction,  since  it  is 
powered  by  the  differential  in  velocity  over  r  [Bourguignon  et  ah,  2013]),  and  a  critical-layer  type  response 
where  locally  Uo(yc)  giving  high  gain.  The  former  effect  is  well  understood  to  lead  to  energy  extraction 
from  the  mean  flow  under  other  formulations.  In  the  latter  case,  the  diagonal  terms  in  the  resolvent  become 
extremely  large  at  the  location  where  c  =  uo(j/)  with  large  Reynolds  number.  The  concept  of  a  critical  layer 
in  wall  turbulence  has  received  much  less  attention,  with  the  notable  exception  of  the  analysis  by  Sreenivasan 
[1988].  This  phenomenon  occurs  when  the  resolvent  at  this  location  becomes  large  because,  at  high  enough 
Reynolds  number,  the  response  at  the  critical  layer  to  forcing  is  almost  singular;  the  relevant  eigenmodes 
approach  neutral  stability  in  the  inviscid  limit. 

A  review  of  the  body  of  work  on  critical  layers  in  laminar  flow  is  given  by  Maslowe  [1986]  and  Schmid  and 
Henningson  [2001].  There  are  some  essential  differences  to  the  current  development:  we  address  the  turbulent 
case,  we  maintain  the  presence  of  the  nonlinear  forcing  term  and  (specific  to  pipe  flow)  the  wall-normal 
velocity  and  vorticity  in  the  Orr-Sommerfeld  formulation  of  the  problem  do  not  decouple  in  cylindrical 
coordinates  [Schmid  and  Henningson,  2001].  However  critical  layer  concepts  are  useful,  and  qualitative 
connections  are  made  between  the  theory  and  turbulent  flow  in  Cartesian  geometries  in  McKeon  and  Sharma 
[2010], 


4.1.4  Closing  the  loop:  an  explicit  treatment  of  the  nonlinearity 


The  linear  selection  mechanism  explains  much  of  the  form  of  the  structure  of  real  flow;  however  the  nonlin¬ 
earity  is  also  important.  Though  it  is  not  the  focus  of  this  work,  we  will  briefly  discuss  it  here  for  the  case 
of  a  discretised  field,  i.e.  discrete  K.  Using  the  resolvent  decomposition  of  equation  7  and  equations  9,  11 
and  12,  expressing  the  fields  in  terms  of  the  coefficients  with  j,  a,  b  indexing  over  all  IK,  we  get 

u  =  E  Xji>j  =  Z  fAi  (13) 

3  3 


where 

fj  =  E  (“Vv  •  VV>b,  (t>j)r  XaXb  =  E  NjabXaXb •  (14) 

a,b  a,b 

Here,  Xj  are  the  coefficients  describing  the  projection  of  the  velocity  field  onto  the  velocity  response  modes, 
Xj  =  <7 j  Xj ,  and  ^fjab  is  the  coupling  between  any  three  traveling  wave  singular  modes.  Note  that  Njab  =  0 
where  Fourier  modes  are  triadically  incompatible.  The  full  nonlinear  problem  can  then  be  reduced  to  the 
solution  of 

Xj  =  E  ®jNjo,bXaXb'  (15) 

a,  b 

The  modeling  question  is  the  magnitude  and  sparseness  of  both  Njab  and  a  j .  Notice  that  high  gain  <jj  can 
compensate  for  low  (but  non-zero)  Njab  and  the  quadratic  property  of  f  can  determine  the  amplitudes.  Since 
(j j  and  Najb  are  calculated  from  the  NSE,  this  system  of  equations  together  with  the  assumed  mean  velocity 
profile  constitutes  a  complete  description  of  turbulence  with  transient  behavior  removed.  Sparseness  makes 
this  representation  very  efficient. 

We  see  that  self-sustaining  mode  combinations  are  possible  if  the  nonlinear  forcing  term  resulting  from  the 
interaction  of  component  velocity  modes  is  not  orthogonal  to  the  forcing  modes  required  to  sustain  it.  We 
could  reasonably  look  for  solutions  to  (15)  to  find  self-sustaining  mode  combinations.  Such  a  truncation 
would  give  a  low-order,  discrete  representation  of  (15).  This  is  the  subject  of  ongoing  work. 

This  view  also  gives  an  interesting  insight  into  the  equivalent  of  the  cascade  in  the  inhomogeneous  spatial 
direction.  In  this  direction,  because  of  the  distribution  across  r  of  the  modes,  there  is  a  non-local  (dis¬ 
tributed)  interaction  with  the  mean  shear.  This  possibility  of  extracting  energy  from  the  velocity  profile  is 
mathematically  captured  by  the  non-normality.  The  critical  layer  response,  by  contrast,  is  localized  in  r. 
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The  triadic  interaction  in  homogeneous  isotropic  turbulence  (or  the  t,  x  and  9  directions)  is  well  understood; 
the  mechanisms  for  energy  transfer  have  now  been  effectively  extended  to  the  inhomogeneous  directions 
by  the  analysis.  The  nonlinearity,  in  addition  to  transferring  energy  across  u>,  k  and  n,  subject  to  the 
rules  of  triadic  interaction,  may  also  transfer  energy  across  response  modes  in  y ,  subject  to  the  three-way 
interaction  Njat,  and  the  singular  values  a j .  The  ability  of  a  mode  to  drive  others  in  the  wall  normal  direction, 
essentially  the  wall-normal  cascade,  is  determined  by  the  localization  of  the  mode  in  that  direction,  which 
in  turn  determines  Njab-  This  concept  is  explored  further  in  Bourguignon  et  al.  [2013]. 


4.1.5  Where  to  look  in  3D  spectral  space? 

In  the  absence  of  a  solution  to  equation  (15)  that  identifies  the  modes  required  to  sustain  the  turbulence 
with  Njab  *  0,  our  analysis  at  this  stage  is  guided  by  the  range  of  K  observed  in  real  flows,  which  has 
been  shown  to  correspond  well  to  the  range  where  the  first  singular  values  have  large  magnitude  by  Moarref 
et  al.  [2013]).  This  localization  and  high  response  at  the  critical  layer  dictates  that  the  turbulent  flow’s 
energy  is  concentrated  around  a  thin  “spine”  in  (k,n,u>,l)  space,  where  l  is  the  singular  value  index,  which 
essentially  describes  a  low-dimensional  attractor  in  this  space.  This  will  also  be  discussed  in  forthcoming 
work.  As  such,  we  briefly  review  work  that  has  characterized  the  scales  and  wall-normal  variation  of  the 
three-dimensional,  spatio-temporal  K  spectrum  (which  in  its  entirety  is  still  a  relatively  uncommon  result) 
and  has  thus  informed  our  work. 

The  statistics  and  spectra  of  the  turbulent  fluctuations  have  been  historically  well  studied,  such  that  results 
including  the  Reynolds-number  independence  of  the  near-wall  cycle  and  the  energetic  importance  of  the  very 
large  scales  emerging  with  increasing  Reynolds  number  are  now  well  established.  However,  the  classification 
of  coherent  structure  and  their  connection  to  velocity  statistics  can  be  said  to  be  a  more  recent  focus, 
significantly  enabled  by  advances  in  computational  and  experimental  diagnostic  technology.  The  essential 
structural  features  of  wall  turbulence  have  been  listed  recently  by  Smits  et  al.  [2011].  In  order  of  increasing 
streamwise  scale,  this  list  consists  of  near-wall  streaks  of  streamwise  velocity,  hairpin  vortices,  large  scale 
motions  (LSM)  believed  to  correspond  to  organization  of  hairpin  vortices  into  packet  structures,  and  very 
large  scale  motions  (VLSM)  with  streamwise  extent  on  the  order  of  ten  times  the  outer  lengthscale.  The 
temporal  frequency  is  most  simply  addressed  through  the  range  of  convection  velocities  observed  in  the  flow 
[Morrison  and  Kronauer,  1969,  Morrison  et  al.,  1971,  del  Alamo  and  Jimenez,  2006,  Dennis  and  Nickels, 
2008,  LeHew  et  al.,  2011]. 

Assuming  the  approximate  universality  of  spectra  between  canonical  flows  documented  by  Monty  and  Chong 
[2009]  and  the  range  of  scales  described  in  Smits  et  al.  [2011]  and  references  therein,  conservative  estimates 
of  the  energetically  active  K  =  ( k,n,c )  ranges  are  then  0.1  <  k  <  2ttR+/100,  27t/0.1  <  n  <  2nR+  /10  and 
10/ug|y+=/j+  <  c  <  1.  Practically,  the  sum  of  K  and  its  complex  conjugate,  and  both  left-  and  right-going 
waves  (±n)  must  be  considered  in  order  to  obtain  real-valued  amplitudes  and  avoid  a  non-zero  mean  spanwise 
velocity,  respectively.  In  what  follows,  we  constrain  our  interrogation  of  the  resolvent  analysis  to  these  values 
and  document  well-known  features  of  wall  turbulence  that  are  captured  by  the  rank-1  approximation  to  the 
resolvent. 

4.2  Reproduction  and  extension  of  statistical  and  structural  results  in  unper¬ 
turbed  wall  turbulence 

An  array  of  quantifiable  features  of  wall  turbulence  have  been  well  described  in  the  extensive  wall  turbulence 
literature,  including  velocity  statistics,  spectral  information  and  organizations  of  coherent  structure,  even 
if  their  Reynolds  number  scaling  and  dynamical  importance  still  remain  sufficiently  elusive  to  preclude 
effective  low  order  modeling.  Of  course,  even  observing  these  features  at  high  Reynolds  number  presents 
severe  experimental  and  numerical  challenges,  but  assimilation  of  the  characteristics  of  all  such  recognizable 
results  is  required  for  a  full  understanding  of  wall  turbulence. 

Our  analysis  naturally  leads  to  information  related  to  velocity  spectra  and  distribution  of  statistics  in  the 
wall-normal  direction,  and  we  will  present  selected  results  in  this  subsection.  We  also  target  the  different 
classes  of  structures  listed  in  subsection  4.1.5  above,  demonstrating  that  the  analysis  has  the  potential  to 
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Figure  18:  (a)  Distribution  of  streainwise  energy  over  the  pipe  radius  for  the  near  wall  mode  with  K  = 
(2irR+/1000,2nR+/100, 10/uJ |y+=jR+ ) .  Reynolds  numbers:  -  75xl03,  —  150xl03,  -  410xl03,  —  lxlO6.  (b) 
Shape  of  the  first  singular  mode  representative  of  the  dominant  near  wall  motions.  Color  denotes  isosurfaces 
of  streamwise  velocity  (streaks) ,  where  red  and  blue  correspond  to  high  and  low  velocity  respectively  relative 
to  the  mean  flow  (heading  into  the  page),  and  the  white  arrows  show  the  sense  of  the  in-plane  velocity  field. 


unify  the  statistical  and  structural  observations  of  wall  turbulence  and  lead  to  inferences  about  the  key 
underlying  physics. 

We  note  that  experimental  studies  typically  interrogate  data  based  on  averaging  conditional  on  some  charac¬ 
teristic  typical  to  the  structure  to  be  examined.  For  instance,  Dennis  and  Nickels  [2011b]  average  conditional 
on  a  swirl  criterion  to  identify  hairpin  vortices,  yielding  the  scale  distribution  in  K  and  y  associated  with 
the  structure.  In  contrast,  we  work  in  the  opposite  direction:  specifying  K  and  yielding  the  wall-normal 
distribution  of  the  velocity  of  the  response  modes,  which  we  associate  with  coherent  structure. 


4.2.1  Near- wall  cycle 

The  structure  and  dominant  dimensions  associated  with  the  autonomous,  near-wall  cycle  of  turbulence  have 
been  known  since  the  seminal  experiments  of  Kline  et  al.  [1967a].  Significant  study  of  the  dynamics  of  the 
flow  in  this  region,  performed  by,  e.g.  Waleffe  [1997],  Jimenez  and  Pinelli  [1999],  Schoppa  and  Hussain  [2002] 
using  low  Reynolds  number  data,  revealed  the  well-known  quasi-streamwise  vortex  and  streamwise  velocity 
streak  structure,  while  experiments  in  the  atmospheric  surface  layer  [Klewicki  et  al.,  1995]  have  demonstrated 
the  Reynolds  number  independence  of  the  streak  spacing  in  viscous  units  at  the  near-wall  scales  identified 
above,  namely  A£  r;  1000  and  A^  rj  100.  Comparison  with  the  structure  arising  from  our  analysis  requires 
the  specification  of  an  additional  parameter,  the  phase  velocity,  c.  For  the  dominant  wavelengths  associated 
with  the  near-wall  cycle,  the  first  singular  mode  is  critical  and  attached  to  the  wall,  as  identified  in  figure  17, 
when  c+  =  10  -  15,  the  minimum  convection  velocity  associated  with  energetic  disturbances. 

The  wall-normal  distribution  of  the  streamwise  component  of  the  first  singular  mode  with 

K  =  (27ri?+/1000, 27ri?+/100, 10/uJ  |y+=jR+ ) 

is  essentially  Reynolds  number-independent,  as  shown  by  the  distributions  of  the  amplitude  over  an  order 
of  magnitude  in  Reynolds  number  in  figure  18(a).  In  fact,  this  independence  emerges  directly  from  the 
Cartesian  coordinate  version  of  critical  layer  scaling  [McKeon  and  Sharrna,  2010]  for  modes  with  K  constant 
in  inner  units  in  a  region  where  the  mean  velocity  profile  is  self-similar  in  inner  units,  i.e.  for  sufficiently  high 
Reynolds  numbers.  The  first  singular  velocity  response  mode  captures  this  scaling  and  thus  postdicts  from 
a  low  rank  approximation  to  the  NSE  that  the  location  of  the  near-wall  cycle  activity  should  be  Reynolds 
number  independent.  Similar  results  can  be  obtained  for  outer  scaling  modes  [McKeon  and  Sharrna,  2010]. 

The  velocity  components  are  distributed  in  such  a  way  as  to  create  a  periodic  array  of  alternating  streamwise 
velocity  streaks  accompanied  by  streamwise-aligned  vortices  in  the  cross-stream  plane,  which  are  elevated 
slightly  from  the  wall  with  increasing  downstream  distance.  See  figure  18(b),  figure  1  from  McKeon  et  al. 
[2010],  for  a  cartoon  of  the  structure  associated  with  this  near- wall  mode.  There  is  a  strong  resemblance 
to  a  system  of  quasi-streamwise  vortices  and  streaks  of  streamwise  velocity,  with  the  simplicity  of  this 
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Figure  19:  Wall-normal  variation  of  (a)  amplitudes  of  the  three  velocity  components  and  (b)  corresponding 
phases  in  multiples  of  2-7T  for  the  VLSM  response  mode  with  K  =  (1, 10,2/3)  at  R+  =  1800.  The  lines  denote 
—o~:  u ,  v,  -  x  w.  Solid/dashed  horizontal  line  shows  the  location  of  the  critical  layer. 


representation  made  possible  because  of  the  absence  of  decorrelation  associated  with  the  many  other  (less- 
energetic)  K  modes  that  are  active  in  this  region  of  the  flow. 

Thus  we  see  that,  first,  the  rank-1  approximation  to  the  resolvent  generates  velocity  response  modes  that 
display  the  scaling,  distribution  and  structure  consistent  with  the  state  of  understanding  concerning  the 
near-wall  cycle.  Secondly,  a  modified  critical  layer  analysis  provides  scaling  consistent  with  observations  of 
wall  turbulence. 


4.2.2  Characteristics  of  the  very  large  scale  motions 

Similar  success  can  be  obtained  by  studying  the  first  singular  velocity  response  mode  at  the  K  combination 
believed  to  correspond  to  the  very  large  scale  motions  (VLSMs).  Figure  19  shows  the  response  mode 
amplitude  and  phase  associated  with  K  =  (1, 10,2/3),  selected  by  McKeon  and  Sharma  [2010]  to  best  match 
a  series  of  observations  of  VLSMs  in  the  literature  that  are  incomplete  in  terms  of  the  three  parameters 
required  for  input  to  our  analysis.  Irrespective  of  the  exact  values  selected,  the  general  form  of  this  mode 
is  similar  to  the  attached  critical  mode  shown  in  figure  16,  namely  streaks  of  streamwise  velocity  with  a 
small  inclination  from  the  wall  and  associated  cross-stream  vortices,  in  agreement  with  the  observations 
of  Hutchins  and  Marusic  [2007]  and  Chung  and  McKeon  [2010].  The  strength  of  each  of  these  features, 
common  to  attached  critical  inodes,  depends  on  the  aspect  ratios  given  by  the  ratios  of  k,  n  and  c,  where 
the  last  dictates  an  approximate  wall-normal  scale  associated  with  the  wall-normal  distance  to  the  critical 
layer.  Consequently,  the  cross-stream  component  is  weak  for  the  VLSMs. 

Of  particular  interest  are  the  distributions  of  the  wall-parallel  components  of  the  mode.  Near  the  wall,  in 
the  region  that  gives  the  mode  its  “attached”  designation,  the  phase  of  the  mode  monotonically  decreases, 
leading  to  the  aforementioned  small  angle  of  inclination  to  the  wall.  However,  further  from  the  wall  the 
amplitude  is  smaller,  but  there  is  essentially  zero  phase  variation  in  the  wall-normal  direction.  Not  only  does 
this  behavior  mirror  the  expected  phase  variation  beyond  the  ir  phase  jump  associated  with  classical  critical 
layer  theory  [Maslowe,  1986],  but  it  is  also  in  good  agreement  with  cross-correlation  results  in  the  literature 
that  had  remained  somewhat  puzzling  [Tutkun  et  ah,  2008,  Guala  et  ah,  2010,  Chung  and  McKeon,  2010]. 
Guala  et  al.  [2010]  observed  similar  shapes  in  the  two-point  cross-correlations,  Rxx{r^.,  y+ ,  y^.ef),  in  the  near- 
thermally-neutral  atmospheric  surface  layer,  with  temporal  records  transformed  to  the  spatial  domain  using 
Taylor’s  hypothesis.  The  expected  near- wall  structural  angle  was  observed  with  a  reference  location,  y^ef, 
close  to  the  wall,  while  a  sufficiently  large  y^ej  led  to  effective  capture  of  a  region  of  weaker  correlation  with 
reduced  phase  variation  far  from  the  wall.  Interestingly,  Chung  and  McKeon  [2010]  found  similar  behavior  of 
the  velocity  Held  in  LES  of  long  channels  conditioned  on  the  occurrence  of  large-scale  negative  and  positive 
velocity  excursions,  with  the  same  dependence  on  the  wall-normal  location  of  the  condition.  Note  that  while 
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Figure  20:  Kelvin’s  cats’  eyes  (also  called  Kelvin-Stuart  vortices):  structure  at  the  inviscid,  laminar  critical 
layer  as  observed  by  an  observer  moving  at  the  critical  velocity. 


there  is  some  ambiguity  concerning  the  interpretation  of  streamwise  scale  from  both  the  cross-correlation 
results  and  conditional  average  procedures,  it  is  clear  that  the  dominant  scales  are  of  the  order  of  several 
boundary  layer  thicknesses  long  in  the  streamwise  direction.  It  also  seems  that  the  general  variation  of  phase 
in  the  wall-normal  direction  described  above  is  a  common  feature  that  can  be  more  easily  isolated  by  these 
techniques  as  the  Reynolds  number  increases  and  the  peak  amplitude  of  the  LSM/VLSM  both  increases  and 
occurs  further  from  the  smaller-scale  near-wall  activity.  We  note  also  that  the  mode  shown  in  figure  19  has 
a  non-negligible  amplitude  over  at  least  30%  of  the  radius,  in  which  distance  the  change  in  mean  velocity  is 
approximately  80%  of  the  centerline  velocity,  underscoring  the  underlying  reason  for  the  failure  of  Taylor’s 
hypothesis  for  very  large  streamwise  scales  [Monty  and  Chong,  2009]. 

For  K  =  (1, 10, 2/3),  a  simple  prediction  of  the  scaling  of  the  peak  location  in  streamwise  amplitude,  y+,  can 
be  made,  where  y +  is  understood  to  track  with  reasonable  fidelity  the  wall-normal  location  of  the  critical 
layer.  McKeon  and  Sharrna  [2010]  showed  that  the  wall-normal  location  where  the  response  mode  with 
(k,n)  =  (1, 10)  first  attains  attached  critical  status  occurs  at  c  =  2/3,  independent  of  Reynolds  number,  and 
that  this  location  lies  within  the  region  of  logarithmic  scaling  of  the  mean  velocity.  The  latter  observation 
allows  expression  of  the  y+  in  terms  of  the  local  mean  velocity  at  the  critical  layer,  two-thirds  of  the  centerline 
velocity,  solely  by  taking  advantage  of  the  log  law  and  Reynolds  similarity  of  the  outer  flow: 

—  \ny+  +  B  =  -  (  —  In R+  +  B  +  c\  .  (16) 

K  3  \K  ) 

Here  k  is  the  von  Karrnan  constant,  and  B  and  C  are  the  additive  constants  associated  with  the  mean 
velocity  in  the  log  region  and  the  constant  wake  function.  Solution  of  equation  16  with  standard  values  for 
the  constants  [McKeon  and  Sharma,  2010]  leads  to  the  prediction  of 

y+  =  0.8i?+  2/3 ,  (17) 

which  appears  to  be  well-borne  out  by  examination  of  the  experimental  data  of  Morrison  et  al.  [2004] 
(McKeon  and  Sharma  [2010],  figure  13). 

Despite  the  difficulties  applying  critical  layer  analysis  to  turbulent  pipe  flow  identified  earlier,  it  is  worth 
mentioning  a  related  result  that  is  at  least  encouraging.  McKeon  and  Sharma  [2010]  observed  that  classical 
critical  layer  analysis  yields  a  scaling  with  R+  2/3  for  the  critical  layer  in  the  upper  branch  solution  for 
neutrally  stable  modes,  along  with  the  more  recognizable  R+  1^2  for  the  wall  layer  in  both  upper  and  lower 
branches.  The  critical  layer  in  the  lower  branch  scales  with  R+  1^5.  That  the  exponent  from  the  classical 
analysis  is  in  agreement  with  the  scaling  of  the  VLSM  energetic  peak  gives  hope  that  the  c  =  2/3  for  the  VLSM 
determined  observationally  from  our  analysis  is  a  fundamental  theoretical  result.  However,  the  comparison 
of  the  variation  of  the  location  of  the  VLSM  peak  in  the  boundary  layer  data  of  Guala  et  al.  [2010]  and  (a 
proxy  for  it)  from  Mathis  et  al.  [2009]  shown  in  figure  16  of  McKeon  and  Sharma  [2010]  with  the  pipe  result 
of  equation  17  gives  a  hint  that  this  behavior  (the  branch  choice)  may  not  be  independent  of  flow  geometry. 

4.2.3  Hairpin  vortices  and  structural  organization 

Our  discussion  of  self-organization  of  coherent  structure  and  its  relationship  with  energetic  large  scale  motions 
afforded  from  resolvent  analysis  begins  with  reference  to  the  classical  structure  of  the  critical  layer  described 
by  Lord  Kelvin,  which  hints  at  the  origins  of  hairpin  vortices.  Figure  20  shows  a  sketch  of  the  periodic, 
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closed  streamline  structure  seen  by  an  observer  moving  at  the  critical  velocity  associated  with  a  classical, 
inviscid  (two-dimensional)  critical  layer.  The  three-dimensional  wavenumbers  expected  from  observations 
of  wall  turbulence  (as  opposed  to  the  two-dimensional  perturbations  that  emerge  as  the  most  unstable  in 
linearly  unstable  laminar  flow  via  Squire’s  theorem)  suggest  that  any  structure  derived  in  the  turbulent 
case  will  be  more  complex:  at  the  very  least,  the  superposition  of  left-  and  right-running  velocity  response 
modes  suggests  a  three-dimensional  vorticity  variation,  which  we  will  show  below  appears  to  represent  the 
commonly-drawn  outline  of  a  hairpin  vortex. 

A  fundamental  concern  with  considering  coherent  vortical  structure  in  a  turbulent  field  associated  with  a 
strong  mean  shear  is  the  challenge  of  objective  identification  of  rotation  over  shear.  Recent  reviews  of  this 
problem  are  given  by,  e.g.,  Jeong  and  Hussain  [1995]  and  Chakraborty  et  al.  [2005];  for  the  simplified  flow 
fields  associated  with  the  linear  superposition  of  velocity  response  modes  investigated  here,  all  common 
identifiers  give  similar  results,  hence  we  choose  to  work  with  swirl,  A,  which  is  given  by  the  magnitude  of  the 
imaginary  part  of  the  complex  conjugate  eigenvalue  pair  associated  with  the  velocity  gradient  tensor.  This 
is,  as  are  all  the  common  diagnostics,  a  nonlinear  function  of  velocity  (a  characteristic  that  leads  rapidly  to 
structural  complexity,  as  we  will  show  below).  Therefore,  while  the  velocity  response  modes  can  be  linearly 
superposed  to  obtain  approximations  to  the  full  velocity  field  of  increasing  veracity,  a  faithful  representation 
of  the  full  swirl  field  can  only  be  obtained  by  summation  of  all  active  velocity  response  modes  (with  the 
correct  amplitudes).  However,  our  aim  here  is  to  show  that  the  recognizable  foundations  of  the  structure 
observed  in  full  flow  can  be  traced  to  resolvent  analysis,  and  specifically  that  phenomena  like  evolving  packets 
of  hairpin  vortices  arise  naturally  from  the  resolvent  analysis  as  studiable  sub-units. 

Figure  21  shows  an  isosurface  of  the  swirl  field  associated  with  K\  =  (6, 6, 2/3),  a  critical  mode  in  the 
terminology  of  subsection  4.1.3  that  is  just  detached  from  the  wall.  Sharma  and  McKeon  [2013a]  selected 
this  response  mode  as  one  of  a  triadically-consistent  set  that  includes  the  VLSM  (we  summarize  the  full 
corresponding  swirl  field  below).  A  periodic  array  of  pro-  and  retro-grade  hairpin-like  vortices  is  found,  where 
prograde  implies  rotation  in  the  classical  sense  with  a  contribution  to  —uv,  and  the  opposite  for  a  retrograde 
vortex.  This  is  as  expected  since  the  response  modes  represent  the  velocity  distribution  relative  to  the  local 
mean  velocity  and  there  can  be  no  mean  contribution  to  the  mean  shear  from  the  response  mode.  These 
coherent  vortical  structures  can  be  understood  phenomenologically  to  be  the  three-dimensional  equivalent 
of  Kelvin’s  cats’  eyes.  Reference  to  figure  17  confirms  that  the  associated  Reynolds  stress  is  localized  around 
the  critical  layer  and  detached  from  the  wall,  making  the  characteristics  of  this  mode  consistent  with  the 
foundational  elements  of  Townsend’s  attached  eddy  hypothesis. 

The  swirl  function  identifies  regions  of  rotational  motion,  but  it  is  susceptible  to  contamination  by  local 
shear  gradients  in  the  sense  that  the  absolute  value  of  the  swirl  can  be  either  increased  or  decreased  by 
the  simultaneous  presence  of  shear.  Most  notably,  the  positive  shear  associated  with  the  mean  profile  can 
be  shown  [Sharma  and  McKeon,  2013a]  to  increase  (suppress)  the  swirl  value  associated  with  prograde 
(retrograde)  hairpin  vortex  heads.  This  is  a  thresholding  problem:  rotational  motion  is  still  captured  for 
both  senses  of  vortex  rotation,  but  the  usual  thresholding  techniques  applied  to  velocity  fields  lead  to  a 
systematic  under-accounting  of  retrograde  vortices.  McKeon  et  al.  [2010]  explored  this  phenomenon  in 
the  context  of  this  work  with  respect  to  observations  of  swirl  in  the  streamwise/wall-normal  plane  of  a  zero 
pressure  gradient  turbulent  boundary  layer,  while  Carlier  and  Stanislas  [2005]  documented  earlier  the  general 
characteristics  of  the  distributions  of  pro-  and  retro-grade  hairpin  vortices. 

In  the  context  of  exploiting  the  resolvent  analysis  to  demonstrate  recognizable  structure,  we  explore  further 
the  effect  of  local  shear  associated  with  combinations  of  velocity  response  modes.  An  energetic  VLSM,  for 
example,  can  be  viewed  as  enforcing  a  periodic,  enhanced  and  reduced  local  Reynolds  number  on  the  smaller 
motions,  accompanied  by  positive  and  negative  shear,  dun/dy.  In  this  sense,  it  can  be  considered  to  be  an 
inactive  motion  in  the  terminology  of  Townsend.  A  similar  effect  occurs  in  the  other  velocity  components 
and  gradient  directions,  but  the  strongest  of  these  is  the  wall-normal  gradient  of  streanrwise  velocity,  at 
least  for  the  large  scales  which  have  long  been  associated  with  vortex  organization  in  the  literature.  Thus 
the  distribution  of  hairpin  vortices  with  a  particular  K,  which  is  periodic  in  the  absence  of  other  modes, 
will  be  modified  by  local  shear  associated  with  the  linear  superposition  of  any  other  (sufficiently  energetic) 
velocity  modes,  and  this  modification  gives  rise  to  the  self-organization  of  vortical  structures  familiar  in  the 
literature. 
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Figure  21:  The  periodic  array  of  pro-  and  retro-grade  hairpin  vortices  associated  with  near- wall  velocity 
response  modes  at  R+  =  1800,  identified  by  an  isosurface  of  constant  swirling  strength  at  50%  of  the  absolute 
maximum  value,  color-coded  with  the  local  (model)  azimuthal  vorticity.  Red  and  blue  denote  rotation  in 
and  counter  to  the  sense  of  the  classical  (prograde)  hairpin  vortex  respectively. 
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This  phenomenon  is  described  at  length  in  Sharma  and  McKeon  [2013a],  but  we  demonstrate  the  complexity 
of  vortical  structure  that  can  arise  under  the  superposition  of  three  modes  in  figure  22.  The  three  modes 
K\  -  (6, 6, 2/3),  I<2  =  (1,6, 2/3)  (VLSM)  and  K3  =  (7, 12, 2/3)  are  triadically-consistent  in  K  and  overlapping 
in  the  y  direction,  and  will  be  phase-locked  in  the  sense  that  the  identical  wavespeed  means  that  there  is  no 
relative  motion  between  the  modes  and  therefore  no  evolution  of  the  packet.  Sharma  and  McKeon  [2013a] 
selected  relative  mode  amplitudes  consistent  with  experimental  observations,  specified  a  7t/2  phase  difference 
between  the  large-scale  envelope  of  the  two  shorter  inodes  and  the  VLSM  in  the  streamwise  sense  and  aligned 
the  peak  streamwise  intensities  in  the  azimuth,  in  order  to  generate  the  complex  swirl  distribution  shown 
in  figure  22,  which  is  recognizably  a  series  of  packets  of  hairpin  vortices.  Three  prograde  vortices  with  two 
retrograde  ones  sitting  downstream  and  above  the  main  packet  can  be  clearly  identified,  along  with  some 
weaker  vortices  of  both  senses  below  the  packet.  The  prograde  hairpin  heads  occur  further  from  the  wall 
with  increasing  streamwise  direction  and  are  collocated  with  the  inclined  interface  between  positive  and 
negative  streamwise  velocity  fluctuation  imposed  by  the  VLSM  (the  location  of  maximum  shear),  giving  the 
packet  its  distinctive  shape.  We  find  the  complexity  of  coherent  structure  arising  from  even  the  idealized 
arrangement  of  only  three  inodes  to  be  striking. 

This  assembly  of  modes  was  deemed  an  “ideal  packet”  by  Sharma  and  McKeon  [2013a]  because  of  its  lack  of 
evolution  (matched  wavespeeds)  and  spatial  alignment.  Note  that  this  is  developed  as  an  example  of  structure 
from  the  resolvent  analysis  that  can  be  simply  interrogated  rather  than  a  structure  that  we  expect  to  observe 
in  a  real  flow.  Relaxing  any  of  the  alignments  or  superposing  more  velocity  response  modes  will  disrupt  the 
coherence  of  such  a  packet,  such  that  isosurfaces  of  the  swirl  diagnostic  identify  the  canes  and  arches  generally 
observed  in  moderate  Reynolds  number  wall  turbulence.  Permitting  a  range  of  wavespeeds  leads  to  packet 
evolution  in  space  and  time.  However,  Sharma  and  McKeon  [2013a]  speculate  that  this  combination  of 
modes  is  dynamically  important  and  likely  to  be  self-sustaining,  in  addition  to  being  triadically  compatible, 
and  termed  it  a  “turbulence  kernel” .  It  should  be  understood  that  the  Reynolds  stress  associated  with  the 
correct  assembly  of  velocity  response  modes  supports  the  (currently  assumed)  mean  velocity  profile,  but  in 
the  resolvent  interpretation  this  comes  from  the  mean  contributions  to  the  —uv  distribution  rather  than  the 
concentration  of  stress  associated  with  the  hairpin  heads,  both  of  which  can  be  observed  in  the  bottom  panes 
of  figure  17.  The  resolvent  analysis  offers  a  different  explanation  for  the  source  of  the  Reynolds  stress  from 
that  offered  by  the  attached  eddy  hypothesis.  In  the  latter,  the  distribution  of  eddies  effectively  describes 
the  mean  velocity  profile  and  the  scaling  of  the  velocity  fluctuations. 

A  simple  superposition  of  first  singular  modes  gives  clear  insight  into  a  mechanism  of  self-organization  in  wall 
turbulence.  In  crude  terms,  the  local  variations  of  shear  associated  with  the  linear  superposition  of  velocity 
modes  bias  the  identification  of  structure  to  the  edges  of  regions  of  energetic  large  scale  streamwise  motion 
(the  zones  of  constant  momentum  of  Adrian  et  al.  [2000]),  implying  that  coherent  vortex  structure  is  most 
likely  to  be  observed  around  streamwise  streaks  of  low  velocity,  as  shown,  e.g.,  in  the  elegant  conditional 
averaging  of  time-resolved  PIV  by  Dennis  and  Nickels  [2011a].  Perhaps,  then  the  well-known  organizational 
characteristics  of  LSMs  and  VLSMs  in  terms  of  the  observed  preferential  wall-normal  and  spanwise  locations 
of  hairpin  vortices  flanking  large  scale  energetic  motion  can  be  at  least  partially  attributed  to  the  choice  of 
nonlinear  diagnostic  for  measuring  swirling  motions. 

4.3  Summary 

The  agreement  between  response  modes  and  observations  of  the  real  flow  detailed  above  suggests  that  the  low 
rank  approximation  to  the  NSE  afforded  by  the  first  singular  response  modes  over  the  range  of  wavenumbers 
and  frequencies  observed  in  real  flows  is  a  useful  tool  for  modeling  wall  turbulence  and  illuminating  dominant 
flow  dynamics.  It  is  important  to  stress  that  the  dynamics  of  the  flow  are  captured  and  explained  by  the 
analysis.  Of  particular  note  is  the  faithful  connection  to  observed  results  and  the  potential  for  informed, 
minimum  effort  control.  The  frequency-domain  interpretation,  essentially  a  decomposition  into  propagating 
waves  in  the  homogeneous  directions,  leads  to  a  self-sustaining  system  in  which  the  Navier-Stokes  nonlinearity 
acts  as  a  feedback  mechanism  to  excite  a  set  of  sub-systems  with  highly  directional  amplification.  The 
formulation  can  appear  somewhat  controversial  because  in  important  respects  it  is  a  linear  analysis  of  a  flow 
that  is  known  to  be  nonlinear,  but,  as  shown  in  figure  13,  this  representation  is  complete  in  the  divergence- 
free  basis.  In  essence,  the  nonlinearity  is  accounted  for  implicitly  through  the  assumption  of  the  turbulent 
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Figure  22:  Isosurfaces  of  constant  swirling  strength  (33%  of  maximum  value)  for  the  “ideal  packet”  at 
R+  -  1800,  color-coded  by  the  local  (model)  azimuthal  vorticity.  Red  and  blue  denote  rotation  in  and 
counter  to  the  sense  of  the  classical  hairpin  vortex  or  prograde  and  retrograde  vortices,  respectively. 


mean  velocity  profile  and  explicitly  as  providing  a  forcing  at  all  wavenumber/frequency  combinations,  which 
is  required  to  sustain  the  turbulent  fluctuations.  More  information  on  the  organization  of  the  flow,  essentially 
the  amplitudes  and  relative  phases  of  the  velocity  response  modes,  appears  to  be  available  from  consideration 
of  the  triadic  interactions  such  as  those  presented  by  our  turbulence  kernel;  this  is  beyond  the  scope  of 
the  current  manuscript,  but  is  explored  in  Sharma  and  McKeon  [2013a].  Another  benefit  of  using  the 
systems  paradigm  is  the  potential  for  carrying  over  closed-loop  control  thinking  and  technology  to  explore 
the  possibilities  for  turbulence  control  in  the  presence  of  an  almost  ideal  model. 

Our  discussion  has  focused  on  the  (linear)  resolvent  operator  relating  the  nonlinear  interaction  between 
scales  and  the  velocity  response  modes  at  each  wavenumber/frequency  combination.  The  resolvent  contains 
an  analytical  description,  derived  from  the  NSE,  of  the  linear  contents  of  Clauser’s  black  box.  To  a  large 
extent,  its  form  determines  the  form  of  the  lower-dimensional  attractor  guiding  turbulent  flow. 

Much  of  the  analytical  development  described  here  has  been  devoted  to  determining  the  appropriate  basis 
for  the  inhomogeneous  wall  normal  direction.  In  the  spirit  of  a  gain  analysis,  we  used  the  singular  value 
decomposition,  which  ranks  input  forcing  distributions  (right  singular  vectors)  at  a  particular  wavenum¬ 
ber/frequency  combination,  I\  ,  in  terms  of  the  amplification  described  by  the  singular  values.  The  observed 
output  will  be  a  product  of  the  amplitude  of  forcing  and  singular  value  summed  over  the  singular  values.  At 
heart,  the  decomposition  yields  a  basis  for  coherence  in  the  wall  normal  direction. 

The  success  of  the  first  singular  velocity  response  modes  (left  singular  vectors)  in  describing  a  range  of 
observed  features  in  unperturbed  high  Reynolds  number  wall  turbulence  demonstrates  that  the  resolvent 
can  be  considered  low  rank  over  the  range  of  K  observed  in  real  flows,  and  rank-1  in  many  cases,  with 
little  error.  In  particular,  the  rank-1  approximation  appears  capable  of  providing  coupled  insight  into  the 
existence  and  organization  of  coherent  structure  within  the  framework  of  the  distribution  of  fluctuating 
velocities,  which  we  believe  to  be  a  first  for  an  NSE-based  analysis. 

We  believe  that  the  resolvent  formulation  represents  a  promising  approach  towards  the  pre-  and  post-diction 
of  the  behavior  of  wall  turbulence,  for  understanding  the  essential  dynamics  of  a  perturbed  flow,  effecting 
control,  and  developing  low  order  models. 
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Limitations  of  the  approach 

We  include  here  a  brief  description  of  the  limitations  of  the  resolvent  analysis,  both  in  its  current  formulation 
and  in  terms  of  extension  to  other  flows. 

A  key  limitation  in  the  current  analysis  lies  in  the  lack  of  information  that  emerges  concerning  relative 
amplitudes  and  phases  between  modes.  This  is  a  consequence  of  the  use  of  unstructured  forcing  and  the 
decomposition  being  optimized  for  “gain” .  Other  choices  of  decomposition  can  give  weight  to  relative  phase 
information,  including  the  sector-bounding  analysis  mentioned  earlierSharma  et  al.  [2011],  Sharma  [2009], 
which  will  be  more  important  where  the  critical  layer  mechanism  is  less  dominant  and  non-normality  drives 
the  modesBourguignon  et  al.  [2013].  This  would  seem  to  correspond  to  regions  where  Taylor’s  hypothesis 
fails.  The  amplitudes  can  be  determined  by  establishing  the  requirements  on  self-sustaining  sub-systems  put 
forward  in  equation  15.  The  velocity  response  modes  contain  information  analogous  to  the  power  spectral 
density,  in  that  amplitude  and  spatial  localization  can  be  determined  by  the  SVD. 

The  first  singular  response  modes  appear  to  capture  many  features  of  wall  turbulence,  but  the  strict  rank-1 
approximation  has  been  shown  to  be  insufficient  to  capture  details  of  the  globally  most  energetic  Fourier 
modes  projected  out  of  DNS  of  turbulent  pipe  flow  [Bourguignon  et  al.,  2013].  These  Fourier  modes  span  a 
significant  portion  of  the  pipe  radius,  leading  to  a  relatively  large  integral  contribution  to  the  turbulent  kinetic 
energy  despite  a  local  amplitude  distribution  that  is  a  fraction  of  that  from  other  K  values.  As  such,  the  wall- 
normal  coherence  of  these  modes  can  be  identified  by  proper  orthogonal  decomposition.  Bourguignon  et  al. 
[2013]  show  that  the  failure  of  the  rank-1  approximation  for  this  class  of  mode  is  associated  with  its  extended 
radial  span  and  can  be  predicted  by  analysis  of  the  resolvent  itself.  As  such,  a  low-rank  approximation  (with 
rank  greater  than  one  but  nonetheless  less  than  ten)  still  appears  to  capture  the  majority  of  the  energy 
content.  We  believe  that  the  absence  of  this  type  of  radially  tall  coherence  in  the  rank-1  approximation  is 
responsible  for  the  absence  of  “bulges”  or  LSM  in  our  analysis.  Again,  we  associate  these  with  modes  for 
which  phase  is  important. 

The  gain  formulation  of  the  SVD  means  that  the  results  are  in  essence  weighted  to  capture  the  most 
energetic  streamwise  velocity  component  for  the  majority  of  K  combinations.  The  success  of  the  rank- 
1  approximation  in  capturing  the  other  velocity  components,  and  in  particular  the  wall  normal  velocity, 
requires  additional  investigation.  The  Reynolds  stress  profiles  associated  with  the  first  singular  response 
show  promising  agreement  with  the  attached  eddy  hypothesis  and  experimental  observations,  but  a  good 
approximation  to  this  distribution  will  be  required  in  order  to  sustain  the  (currently  assumed)  mean  velocity 
profile. 

With  respect  to  application  of  our  analysis  to  other  canonical  flows,  the  spatial  inhomogeneity  associated 
with  a  boundary  layer  also  suggests  the  possibility  of  fundamentally  different  behavior  to  fully-developed 
flows,  particularly  if  the  locally-parallel  approximation  is  not  good. 

Lastly,  we  note  that  the  current  level  of  sophistication  in  this  approach  does  not  support  modeling  based  on 
coherent  structure,  a  method  that  has  been  a  focus  of  research  aimed  at  finding  low  order  models  for  wall 
turbulence  motivated  by  the  importance  of  structure  to  momentum  transport,  concentration  of  Reynolds 
stress,  etc.  Our  work  to  date  identifies  the  complexity  of  coherent  structure,  associated  with  a  nonlinear 
diagnostic,  in  an  otherwise  linear  model.  In  essence,  a  full  description  of  all  energetically-active  modes  may 
be  required  to  observe  the  full  dynamics  and  evolution  of  coherent  structure,  which  would  probably  preclude 
simple  low  order  modeling. 


Future  trends 

We  conclude  by  briefly  outlining  possible  future  developments  ranging  from  further  interrogation  and  mod¬ 
eling  of  the  current  formulation  to  extension  to  account  for  increasing  flow  complexity.  We  believe  that 
features  of  wall  turbulence  beyond  those  described  in  subsection  4.2  are  amenable  to  description  using  our 
analysis.  For  example,  there  appears  to  be  a  strong  relationship  between  the  wall-normal  profiles  of  skewness 
and  the  amplitude  modulation  coefficient,  in  particular  with  respect  to  triadically-consistent  K  combinations 
that  include  the  VLSM  motion  [Sharma  and  McKeon,  2013a].  These  two  measures  can  then  be  interpreted 
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as  giving  information  about  the  relative  phases  between  velocity  response  modes,  one  of  the  missing  pieces 
of  our  analysis  to  date,  as  identified  above. 

We  have  also  noted  the  likely  potential  of  the  resolvent  analysis  to  treat  rough  walls,  or  at  least  those 
with  periodic  geometries.  Under  this  scenario,  admitting  a  spatial  variation  of  the  mean  flow  and  perhaps 
excitation  of  a  range  of  frequencies  due  to  local  flow  separation  seem  to  be  conceptually  simple  extensions 
likely  to  give  insight  into  roughness-induced  changes  to  the  flow. 

An  important  next  step  will  be  to  “close  the  loop” ,  in  the  sense  of  determining  the  amplitudes  of  the  velocity 
response  modes  observed  in  (projected  out  of)  real  flows  that  are  required  to  fulfill  the  condition  of  self¬ 
sustenance  implied  by  the  full  NSE  formulation  of  figure  13.  Once  these  amplitudes  have  been  identified,  a 
full  description  of  energy  transfer  pathways  and  scaling  with  Reynolds  number  can  be  obtained.  Extraction 
of  energy  from  the  mean  flow,  spatial  transport,  spectral  transfer  and  local  triadic  interactions  responsible 
for  sustaining  the  flow  all  follow  from  examination  of  the  resolvent  and  energy  gradients.  The  existence 
of  self-sustaining  turbulence  kernels,  or  limited  combinations  of  modes  that  are  capable  of  self-exciting  via 
triadic  interaction,  is  a  topic  of  current  study.  Well-known  examples  exist  in  the  literature  without  the 
coherent  propagating  wave  interpretation,  but  additional  such  solutions  seem  likely  [Sharma  and  McKeon, 
2013a].  We  note  also  the  likely  connection  between  this  kind  of  development  and  nonlinear  traveling  wave 
solutions  identified  for  transitional  pipe  flow  [Kerswell,  2005].  Closed- loop  control,  given  a  sufficiently  good 
model,  also  seems  a  logical  next  step. 

As  noted  in  the  previous  subsection,  our  analysis  begins  by  assuming  a  mean  velocity  profile,  which  must  be 
sustained  by  the  combined  action  of  the  perturbations  about  it.  Although  the  resolvent  associated  with  the 
mean  velocity  is  contained  in  the  formulation  via  the  Kq  mode,  a  means  of  approximating  it  and  identifying 
the  sensitivity  of  the  response  modes  to  deviation  from  the  true  profile  would  be  useful,  to  say  the  least. 
Towards  this  end,  one  of  the  authors  has  investigated  a  two-dimensional,  three-component  (“2D/3C”)  model 
which  appears  to  isolate  the  appropriate  dominant  nonlinear  mechanism  governing  the  mean  velocity  [Gayme 
et  al.,  2010,  Bourguignon  and  McKeon,  2011]  Future  advances  seem  feasible.  In  this  vein,  we  note  that  all  the 
turbulent  mean  profiles  we  have  encountered  lead  to  stable  Ck ■  With  appropriate  tribute  to  Malkus  [1956], 
we  conjecture  that  the  turbulent  mean  velocity  profile  Uo  always  leads  to  stable  Ck,  to  avoid  a  singular 
resolvent  response. 
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5  Experimental  measurements  of  the  frequency/wavenumber  spec¬ 
trum  and  coherent  structure  using  time-resolved  PIV 

One  component  of  the  work  performed  under  this  grant  was  a  first  characterization  of  the  three-dimensional 
spectrum  of  the  velocity  fluctuations  in  a  canonical  turbulent  boundary  layer,  i.e.  the  decomposition  of  the 
fluctuation  energy  by  streamwise  and  spanwise  wavenumber  and  temporal  frequency.  Typical  experiments 
and  simulations  consider  either  the  streamwise  wavenumber  or  temporal  frequency  spectrum;  conversion 
between  the  two  can  be  estimated  using  Taylor’s  hypothesis  of  frozen  turbulence.  While  some  attention 
has  been  paid  to  the  associated  implications  for  structure  convection  velocity,  there  has  been  only  limited 
investigation  into  the  distribution  of  energy  in  this  spectral  space.  By  employing  time-resolved  particle 
image  velocimetry  (PIV)  in  a  water  channel,  the  measurements  reported  here  detail  the  three-dimensional 
spectrum  in  three  wall-parallel  planes  at  wall-normal  heights  corresponding  to  different  scaling  regions  of  the 
boundary  layer.  Besides  providing  a  fundamental  contribution  to  characterization  of  wall  turbulence,  this 
information  provides  direct  input  to  the  resolvent  analysis  described  in  subsequent  sections. 

In  addition,  the  time-resolved  information  was  used  to  identify  and  track  individual  coherent  structures  in 
the  measurement  domains  and  characterize  their  trajectory,  lifetime  and  footprint  at  the  wall. 

Full  details  of  this  work  reported  in  this  section  were  published  in  Experiments  in  Fluids  as  LeHew  et  al. 
[2011]  and  LeHew  et  al.  [2013]. 

5.1  Experimental  Setup 

Measurements  of  a  zero-pressure  gradient  boundary  layer  developing  over  a  flat  plate  were  performed  on  a 
1.1m  long  by  0.45m  wide  Plexiglass  plate  with  an  elliptical  leading  edge  and  an  adjustable  wedge-shaped 
trailing  edge  flap  situated  in  a  2m  long  by  lrn  wide  test  subsection  of  a  free  surface  water  tunnel  facility. 
The  13.9nnn  thick  boundary  layer  on  the  bottom  of  the  plate  was  studied  to  avoid  interaction  with  surface 
waves.  A  shroud  with  a  wedge-shaped  leading  edge  was  placed  along  the  sides  of  the  plate  to  promote  a 
two-dimensional  flow.  A  thin  strip  of  tape  (approximately  3nnn  long,  0.25mnr  thick,  and  spanning  the  width 
of  the  plate)  was  placed  immediately  downstream  of  the  elliptical  leading  edge  on  the  measurement  side  to 
promote  transition  and  favor  the  establishment  of  a  turbulent  regime  at  the  measurement  location.  The  flow 
was  conditioned  by  passing  it  through  a  perforated  plate,  a  honey  comb,  three  turbulence  reducing  screens 
and  finally  a  6:1  contraction.  Measurements  were  made  starting  0.63m  downstream  of  the  leading  edge. 
The  free  stream  velocity  was  0.67m/s  and  the  free  stream  turbulence  level  was  less  than  0.1%.  The  free 
stream  velocity,  and  thus  Reynolds  number,  for  these  experiments  was  set  to  provide  the  best  flow  conditions 
possible  and  to  allow  a  small  enough  displacement  per  frame  for  accurate  PIV  mesurements.  It  should  be 
noted  that  this  choice  led  to  a  relatively  low  Reynolds  number,  ReT  =  470,  which  limits  the  scale  separation 
observable  in  these  experiments.  A  diagram  of  the  tunnel  and  a  photograph  showing  the  test  subsection  and 
coordinate  system  are  presented  in  Fig.  23. 

Time  resolved  2D  PIV  was  used  to  measure  the  flow  field  in  the  wall-normal  and  wall-parallel  directions  in 
two  separate  experiments  using  a  LaVision  PIV  system  with  2  Photron  Fastcam  APX-RS  high  speed  cameras 
equipped  with  Tamron  SP  AF  180mm  F/5.5  macro  lenses.  The  cameras  acquired  images  at  2000  frames 
per  second  with  1024x1024  pixel  resolution  for  all  experiments.  The  flow  was  seeded  with  10/irn  hollow 
glass  spheres  with  a  specific  gravity  of  1.1.  The  seeding  density  was  similar  in  both  wall  normal  and  wall 
parallel  measurements  and  the  average  particle  image  diameter  was  2.2  and  3.1  pixels,  respectively.  The  flow 
was  illuminated  by  a  Photonics  DM20-527  solid  state  laser  providing  20mJ /pulse  with  a  sheet  thickness  of 
approximately  lnnn.  The  sheet  was  thick  enough  to  nearly  eliminate  out-of-plane  loss  with  \v\At/yo  <  0.04, 
much  less  than  the  maximum  allowable  value  of  0.25  suggested  by  Keane  and  Adrian  [1990],  where  yo  is  the 
laser  sheet  thickness,  At  is  the  time  between  images,  and  v  is  the  out-of-plane  velocity  estimated  from  the 
rms  velocity  fluctuations.  The  two-dimensionality  of  the  flow  in  wall  parallel  planes  is  verified  in  Fig.  24 
showing  the  velocity  averaged  over  the  streamwise  direction  and  time  for  y+  =  34.  For  this  plane,  the  local 
mean  streamwise  velocity  at  a  given  spanwise  location  varied  from  the  global  mean  by  no  more  than  1%. 
The  spanwise  variation  of  velocity  is  similar  at  the  other  wall  normal  locations. 
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(b)  Test  subsection 

Figure  23:  (a):  Water  tunnel  schematic  showing  the  location  of  the  test  subsection,  (b):  Photograph  of 
the  test  subsection;  red  lines  outline  the  plate  and  blue  lines  outline  the  shroud  and  its  wedge-shaped  leading 
edge.  A  top-down  view  is  provided  in  the  inset  to  clarify  the  shroud  orientation. 
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Table  1:  Experimental  parameters,  spectral  measurements 


ReT 

Reg 

S  9  Uoo 

UT 

v/uT 

A  y+ 

A  t+ 

470 

1280 

13.9mm  1.72mm  0.67m/s 

0.03m/s 

29.5  fim 

10.0 

0.5 

Plane 

y+ 

y/S 

U+ 

U/Uoe 

u+ 

^ rms 

A;c+,  Az+ 

Lx/S 

Lz/6 

1 

34 

0.07 

13.1 

0.59 

2.10 

37.6 

9.98 

4.91 

2 

108 

0.23 

16.4 

0.73 

1.53 

37.1 

9.85 

4.84 

3 

278 

0.59 

20.2 

0.90 

1.25 

37.5 

10.03 

4.89 

5.1.1  Measurements 

Wall  normal  measurements  were  taken  first  to  characterize  the  flow  at  the  start  of  the  measurement  location 
using  only  one  camera.  The  field  of  view  was  50.6mmx50.6m?n  (LxxLy).  Due  to  camera  memory  limitations, 
a  maximum  of  2048  images  could  be  taken  at  a  time.  For  collapse  of  the  mean  profile,  5  experiments  were 
performed. 

For  wall  parallel  measurements,  two  cameras  were  placed  side-by-side  in  the  streamwise  direction  with  a 
10-20  pixel  overlap  producing  a  total  field  of  view  of  approximately  140mm  x  70 mm  ( Lx  xLz)  with  some 
slight  deviation  between  the  three  planes  investigated.  A  laser  sheet  was  guided  into  the  test  subsection 
parallel  to  the  flat  plate  and  was  centered  at  locations  1mm,  3.2mm,  and  8.2mm  from  the  wall  in  order  to 
capture  spectral  information  in  the  near  wall,  outer  log  layer,  and  wake  regions  of  the  turbulent  boundary 
layer.  For  each  plane,  40  experiments  were  performed  providing  over  80,000  instantaneous  realizations 
allowing  collapse  of  the  calculated  power  spectra.  Experimental  conditions  for  this  set  of  experiments  are 
summarized  in  Table  1. 

A  second  set  of  wall-parallel  measurements  were  taken  at  four  wall-normal  locations  independently  ranging 
from  y/S  =  0.08-0.48  in  order  to  track  coherent  structures.  Two  cameras  were  placed  side-by-side  in  the 
streamwise  direction  to  elongate  the  streamwise  field  of  view.  25  runs  of  2048  images  were  recorded  at 
2500  Hz,  providing  950  eddy  turnover  times  (Te)  of  data  per  point.  Additional  experimental  parameters  are 
provided  in  Table  2  where  Aa’  and  Az  are  the  vector  spacings  in  the  x  and  z  directions,  respectively,  Lx  and 
Lz  are  the  size  of  the  field  of  view  in  the  x  and  z  directions,  respectively,  At  is  the  time  between  images, 
and  Tm  is  the  total  measurement  time. 

Given  that  most  of  the  swirling  coherent  structures  identified  in  this  study  did  not  extend  more  than  75-90d„ 
(5-6Aa;)  in  either  x  or  z,  the  resolution  of  the  equivalent  radius,  req  and  inclination  angle,  a,  was  limited,  as 
will  be  discussed  in  Subsection  5.4.2.  In  addition,  the  field  of  view,  Lx/5  «  4.3  limited  the  average  observation 
time  of  a  vortex,  T0  =  Lx/U(y)  to  5-6. 3Xe,  with  T0  increasing  with  decreasing  y,  as  will  be  discussed  in 
Subsection  5.5.2. 

Table  2:  Experimental  parameters  for  swirling  coherent  structure  experiments 


ReT 

Reg 

<5  (mm) 

0(mnr) 

Uoc  (m/s) 

uT  (nr/s) 

410 

1085 

13.5 

1.60 

0.63 

0.028 

32.9 

Plane 

y+ 

y/s 

U/Uoe 

Ax+,  A z+ 

Lx/S 

Lz/S 

A  t+ 

Tm/Te 

1 

33 

0.08 

0.63 

15 

4.3 

2.2 

0.34 

950 

2 

58 

0.14 

0.69 

15 

4.3 

2.2 

0.34 

950 

3 

122 

0.30 

0.77 

15 

4.4 

2.2 

0.34 

950 

4 

198 

0.48 

0.86 

14 

4.3 

2.1 

0.34 

950 
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Figure  24:  The  spanwise  variation  of  the  velocity  averaged  over  the  streamwise  direction  and  time  at 
y+  =  34 

5.1.2  Vector  Processing 

The  velocity  field  calculations  for  both  the  wall  normal  and  wall  parallel  planes  were  performed  using 
LaVision’s  Davis  software.  In  the  case  of  the  wall  parallel  planes,  the  vector  fields  were  stitched  together 
after  processing.  To  assure  an  exact  overlap  of  the  processed  fields,  the  original  images  were  cropped  before 
processing  to  ensure  overlap  of  an  integer  number  of  interrogation  windows  after  processing.  In  addition, 
due  to  bubbles  and  particles  occasionally  aggregating  at  the  wall,  for  the  plane  nearest  the  wall,  an  average 
over  all  images  was  taken  and  subtracted  from  each  raw  image  before  vector  calculation  to  reduce  the  effects 
of  stationary  tracers  on  this  calculation. 

A  50%  overlap  of  interrogation  windows  was  used  to  satisfy  the  Nyquist  criterion,  which,  as  will  be  discussed 
in  Sect.  5.3.1,  is  necessary  to  resolve  the  spectrum  properly.  The  correlation  peak  fitting  was  done  using 
a  Gaussian  three-point  estimator  in  each  coordinate  direction  independently,  Whittaker  reconstruction  was 
used  for  image  reconstruction,  and  in-plane  particle  pair  loss  was  minimized  by  using  a  window  shift.  For  wall 
normal  measurements  the  window  size  was  reduced  from  16x16  to  12x12  pixels  over  two  passes  and  for  wall 
parallel  measurements,  both  passes  were  performed  with  a  32x32  pixel  window  with  each  camera  processed 
separately.  The  non-dimensional  vector  spacing,  field  size,  and  wall  normal  locations  are  summarized  in 
Tables  1  and  2. 

After  the  initial  vector  calculation,  spurious  vectors  were  removed  and  replaced  via  interpolation.  In  all 
images,  less  than  3%  of  the  vectors  were  removed  so  the  effect  of  interpolation  on  the  measurements  was 
considered  minimal.  While  the  wall  normal  data  was  smoothed  with  a  3x3  filter,  no  smoothing  was  performed 
in  wall  parallel  planes  as  to  preserve  the  fluctuating  velocity  signal  and  spectral  shape. 

5.2  Flow  Properties  and  Statistics 

The  mean  velocity  profile  and  boundary  layer  thickness  were  calculated  from  the  wall  normal  measurements 
averaged  in  the  nearly  homogeneous  streamwise  direction,  over  all  frames,  and  over  all  experiments.  The 
boundary  layer  thickness,  calculated  using  the  U ( y )  =  0.99t/oo  criterion,  varied  by  5%  over  the  50mm  window. 
The  mean  profile  is  compared  to  the  experiments  of  DeGraaff  and  Eaton  [2000]  at  a  similar  Reynolds  number 
( Reg  =  dUoo/u  =  1430)  in  Fig.  25a.  6  is  the  momentum  thickness  and  Uoo  is  the  free  stream  velocity. 

Since  no  direct  measurements  of  the  wall  shear  stress  were  made  in  this  experiment,  the  friction  velocity  was 
calculated  from  a  least  squares  fit  to  the  predetermined  log  layer  marked  with  circles  in  Fig.  25a. 

In  the  wall  parallel  measurements  S  changed  by  about  10%  over  the  whole  field  of  view.  While  the  outer 
scale,  y/5,  is  not  constant  over  the  measurement  plane,  the  inner  scaled  wall  distance  is  nearly  constant  since 
uT  only  decreases  weakly  with  streamwise  extent.  Thus,  each  plane  of  the  wall  parallel  data  is  representative 
of  a  particular  inner-scaled  distance  from  the  wall.  Data  will  also  be  presented  using  an  averaged  outer  scale 
in  some  cases  for  convenience. 
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(a)  Mean  profile  (b)  Streamwise  rms  fluctuations 

Figure  25:  In  (a):  the  mean  velocity  profile  is  presented  where  the  symbols  are  — Mean  profile  from  wall 

normal  measurements,  ■:  mean  velocity  from  wall  parallel  plane  measurements, - :  Spalding  fit  [Spalding, 

1961],  data  from  DeGraaff  and  Eaton  [2000]  at  Reg  =  1430,  and  O:  points  selected  to  represent  the 
log  layer.  In  (b):  black  lines  and  symbols  are  for  uxms  and  blue  lines  and  symbols  are  for  v+ms  where,  ■: 
u1rms  from  wall  parallel  plane  measurements,  -  •  -  :  v+ms  from  the  wall  normal  plane  measurement,  —  0  — : 
data  from  DeGraaff  and  Eaton  [2000]  at  Reg  =  1430,  and  —  A  — :  data  from  Erm  and  Joubert  [1991]  at 
Reg  = 1003 


The  rms  streamwise  velocity  fluctuations  for  each  wall  parallel  plane  are  presented  in  Fig.  25b  and  match 
closely  with  the  data  from  DeGraaff  and  Eaton  [2000]  and  Erm  and  Joubert  [1991].  Any  small  deviation 
from  the  literature  values  may  be  from  the  spatial  averaging  present  in  the  PIV  calculations  as  discussed 
in  Sect.  5.3.1.  The  wall  normal  fluctuations  measured  in  the  wall  normal  plane  are  also  presented  for 
completeness  and  deviate  from  the  data  shown,  particularly  near  the  wall.  This  deviation  arises  from  limited 
resolution  in  the  wall  normal  plane  and  is  also  observed  in  the  streamwise  velocity  fluctuations  in  this  plane. 
As  already  noted  this  is  not  an  issue  in  the  wall  parallel  planes  which  will  be  the  focus  of  all  discussion  from 
here  on. 


5.3  Spectral  energy  distributions 

5.3.1  Calculation  of  the  power  spectral  density  in  wall  parallel  planes 

The  power  spectral  density  is  defined  as  the  Fourier  transform  of  the  correlation  function  of  the  velocity 
fluctuations,  where  the  relevant  non-normalized  3D  correlation  function  is  defined  in  equation  18. 


RuiUj(Px,Pz,T )  =  (u(x,z,t)u(x  +  px,z+  pz,t  +  r))  (18) 

In  the  above  equation,  px  and  pz  are  the  separation  between  points  in  the  streamwise  and  spanwise  directions, 
respectively,  r  is  a  separation  in  time,  the  subscripts  i  and  j  represent  velocity  components,  and  the  triangle 
brackets  represent  an  ensemble  average.  The  Fourier  transform  pair  relating  the  correlation  function  to  the 
power  spectrum  0  is  given  in  Equations  19  and  20. 


oo 

QuiUj(kx,kz,cv)  =  -^- 3  JJJ  RUtuAP^Pz’T)e-^k*+p*k*+T^dpxdpzdT 

—  oo 

oo 

RutUj(px,Pz,T)  =  fff  Q^U]{kx,kz,u:)e^k^p‘k^Ukxdkzd^ 


(19) 


(20) 
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kx  and  kz  are  the  streamwise  and  spanwise  wavenumbers,  respectively,  and  u>  is  the  angular  frequency. 
The  Fourier  transform  pair  is  defined  as  above  to  allow  for  proper  normalization  of  the  spectrum,  where 
the  normalization  is  defined  by  combining  Equations  18  and  20  at  zero  shift  in  time  and  space  as  given  in 
Equation  21. 


Ru 


oo 

,  (0,0,0)  =  (uiUj)  =  JJJ  (-)UtUj  (kX)kz,oj)dkxdkzd/jj 


(21) 


In  practice,  the  auto-spectrum  is  calculated  as  the  square  magnitude  of  the  finite  Fourier  transform  of  the 
velocity  fluctuations  and  is  normalized  using  Equation  21. 

2D  spectra  can  be  calculated  by  integration  of  the  3D  spectrum  over  one  dimension  as  shown  in  Equation 
22  (where  the  limits  are  finite  in  practice),  or  by  calculating  the  finite  Fourier  Transform  of  a  2D  subset  of 
the  original  3D  sample  records  and  taking  an  ensemble  average  over  all  subsets. 


®UiUj(ki,k2)  =  /  QUiuj(ki,k2,k3)dk3  (22) 

J  —  OO 

k\,  k‘2 ,  and  k3  are  any  three  wave  vectors.  Similarly,  ID  spectra  can  be  calculated  as  the  integral  of  a  2D 
spectrum  or  by  using  a  ID  subset  of  the  original  3D  dataset  and  ensemble  averaging  over  all  subsets. 

As  discussed  by  Adrian  [1988],  the  velocity  measured  by  the  PIV  algorithm,  Ui(x,y,z,t),  is  not  the  velocity  at 
a  single  point  in  the  flow,  Ui(x ,  y,  z,  t ),  but  the  velocity  signal  averaged  over  the  laser  sheet  thickness  (~lmm) 
and  convolved  with  a  rectangular  window,  h(x,z ),  which  represents  the  averaging  effect  of  an  interrogation 
window.  The  relation  between  the  measured  and  the  true  velocity  is  given  in  Equation  23  and  the  measured 
fluctuating  velocity  is  defined  in  Equation  24. 


Ui(x,y,z,t )  = 


X 


y+yo/z 


y~yo/z 


Ui(x,y',  z,t)dy' 


h(x,  z) 


(23) 


Ui(x,y,z,t )  =  Ui(x,y,z,t )  -  Ui(y)  (24) 

U  is  the  instantaneous  velocity,  u  is  the  fluctuating  velocity,  the  star  represents  convolution,  and  the  overbar 
indicates  an  average  over  all  space,  time,  and  experiments  that  closely  approximates  the  true  mean  as 
discussed  by  Tomkins  and  Adrian  [2005]. 

The  spatial  averaging  from  the  interrogation  window,  h,  greatly  attenuates  fluctuations  smaller  than  two 
window  widths,  and  will  also  attenuate  scales  near  this  limit  to  a  lesser  extent.  This,  then,  leads  to  an 
attenuation  of  the  measured  mean  square  fluctuations  and  the  velocity  spectrum. 

When  considering  the  measured  velocity  signal,  u ,  given  in  Equation  24,  a  different  spectrum  will  result.  The 
square  magnitude  of  the  Fourier  transform  of  Equation  24  is  given  in  Equation  25,  where  the  convolution 
theorem  has  been  used  to  separate  out  the  effects  of  smoothing  from  the  interrogation  window. 


kz,ix)  —  (kx,  kz ,  ca)r(/ca,,  /c~) 

r(kx,kz)  =  sine2  ^^kx^  sine2  ^^k^  '  ^ 

0  is  the  3D  power  spectrum  of  the  measured  velocity  signal,  T  is  the  square  magnitude  of  the  Fourier 
transform  of  the  rectangular  window,  h,  and  W  is  the  PIV  interrogation  window  size  in  either  x  or  z  as 
denoted  by  the  subscript.  The  measured  spectrum,  0,  should  be  normalized  by  the  measured  mean  square 
value  of  the  velocity  fluctuations,  which  is  the  mean  square  value  of  the  fluctuations  defined  in  Equation  24. 

The  effect  of  T  is  to  attenuate  the  spectrum,  particularly  at  high  wavenumbers,  where  the  amplitude  of  T 
goes  to  zero  at  k  =  2ir/W.  Since  T  is  known,  one  could  recover  the  original  spectrum,  0  simply  by  dividing  by 
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r,  although,  in  practice,  this  does  not  work  where  T  approaches  zero.  As  noted  by  Foucaut  et  al.  [2004],  the 
attenuation  become  significant  (more  than  50%  attenuation)  above  a  wavenumber  kcut  =  2.8 /W.  Following 
their  criterion,  all  data  beyond  this  cutoff  will  be  ignored. 

Since  the  experimental  data  is  not  periodic  and  of  finite  length,  the  data  is  windowed  prior  to  calculating 
the  Fourier  transform  to  prevent  spectral  leakage.  A  3D  Hanning  window,  an  extension  of  the  2D  Hanning 
window  used  by  Tomkins  and  Adrian  [2005],  is  used  here.  With  the  current  data,  the  effect  of  leakage  is  most 
prevalent  at  the  low  frequency  end  of  the  spectrum.  Multiplying  the  measured  signal  by  the  3D  Hanning 
window,  denoted  by  g,  provides  an  estimate  of  the  spectrum  as  presented  in  Equation  26. 

QuiUj(kx,kz,oj)  =  [[(UiH)  *  G]  [( U*H *)  *  G*]  ]  (26) 

The  hat  indicates  that  this  is  an  estimate  of  the  spectrum  since  the  signal  is  windowed  and  of  finite  length, 
the  capital  letters  are  the  finite  Fourier  transform  of  the  corresponding  lower  case  quantity  where  UU*  =  0, 
HH*  -  F,  and  the  asterisk  represents  the  complex  conjugate.  Estimates  of  the  2D  and  ID  spectra  can  be 
calculated  as  before  and  the  mean  square  fluctuations  are  the  same  as  for  the  non-windowed  data  using  a 
proper  weighting  of  the  windowing  function. 

The  division  of  this  estimate  of  the  spectrum  by  T  to  correct  the  high  wavenumber  range  gives  a  reasonable 
estimate  of  the  true  spectrum  since  the  convolution  with  G  in  spectral  space  mainly  affects  the  low  frequen¬ 
cy/wavenumber  range.  The  same  operations  were  performed  by  Tomkins  and  Adrian  [2005].  In  presenting 
the  spectra,  both  the  corrected  form  (divided  by  T)  and  the  uncorrected  form  will  be  shown  to  illuminate 
the  effects  of  this  correction. 

To  avoid  spatial  aliasing,  interrogation  windows  with  50%  overlap  were  used  which  attenuate  aliased  energy 
content  above  the  Nyquist  frequency.  To  avoid  temporal  aliasing,  the  data  was  oversampled  during  recording, 
low  pass  filtered,  and  then  subsampled  before  calculating  the  spectrum  as  discussed  by  Wernet  [2007].  The 
low  pass  filter  cutoff  for  the  time  signal  is  set  so  that  cocut  =  U{y)kcut. 

To  summarize,  small  scale  features  of  the  flow  are  attenuated  by  the  averaging  over  the  PIV  interrogation 
window  and  large  scale  features  are  smoothed  out  by  spectral  leakage.  To  account  for  these  issues  and 
prevent  aliasing,  the  data  is  recorded  at  a  higher  sampling  rate  than  necessary  and  the  vector  fields  are 
calculated  using  50%  overlap  of  the  interrogation  windows.  Before  any  calculations  are  performed,  the  data 
is  low  pass  filtered  in  time  to  get  rid  of  unwanted  noise.  Prior  to  calculation  of  the  spectrum,  the  signal  is 
split  into  subsets  (split  twice  in  z  and  thrice  in  time  with  each  region  overlapping  by  50%)  and  windowed 
using  a  Hanning  window  to  help  reduce  spectral  leakage.  Finally,  the  finite  Fourier  transform  of  each  subset 
is  calculated,  all  calculations  are  averaged,  and  the  data  is  divided  by  the  smoothing  function  T  to  get  an 
estimate  of  the  true  spectrum  0  where  data  beyond  the  cutoff  in  kx  and  kz  is  ignored.  ID  and  2D  spectra 
are  calculated  by  integration  of  the  resulting  3D  spectrum. 

5.3.2  Spectral  symmetries 

We  present  here  a  discussion  of  the  symmetries  present  in  the  ID,  2D,  and  3D  streamwise  velocity  spectra 
and  the  proper  normalization  of  each  before  moving  on  to  present  the  ID,  2D,  and  3D  spectra  calculated  from 
the  PIV  measurements.  Due  to  a  lower  than  desirable  dynamic  range  of  the  spanwise  velocity  fluctuations, 
neither  the  spanwise  velocity  spectrum  nor  the  streamwise-spanwise  velocity  cross-spectrum  will  be  presented 
here.  For  this  reason,  in  the  remainder  of  this  document,  the  UiUj  subscript  for  the  spectra  will  be  dropped 
since  all  spectra  will  be  streamwise  velocity  spectra. 

The  3D  auto-spectrum  is  even  and  real  which  provides  symmetry  about  any  coordinate  plane  if  the  spectrum 
is  rotated  by  180°  about  the  origin,  and  thus  only  4  independent  octants  exist.  Similarly,  all  2D  spectra 
are  even  and  real  with  rotational  symmetry  about  the  origin,  and  thus  have  only  2  independent  quadrants 
each.  Therefore,  any  2D  spectra  can  be  represented  by  any  half  plane  and  the  3D  spectrum  by  any  4 
connected  octants.  In  this  half  space  representation,  the  magnitude  of  the  spectra  are  doubled  to  conserve 
energy.  The  2  non-redundant  quadrants  of  each  2D  spectrum  are  shown  in  Figs.  26a,  26b,  and  26c  and  the 
four  non-redundant  octants  (choosing  w  >  0)  of  the  3D  spectrum  are  shown  in  Fig.  27  for  one  wall  normal 
location. 
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Figure  26:  The  panels  show  the  half  space  representation  of  each  2D  spectrum  at  y+  =  108.  The  spectra  at 
other  wall  normal  locations  are  qualitatively  similar.  The  levels  shown  are  from  20%  to  80%  of  the  maximum 
energy  in  10%  increments  moving  from  light  to  dark  shades 
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Figure  27:  The  half  space  representation  of  0 (kx,kz,oj)  at  y+  =  108  is  shown,  where  the  spectra  in  other 
planes  are  qualitatively  similar.  The  levels  are  from  15%  to  75%  of  the  maximum  energy  in  20%  increments 
moving  from  light  to  dark  shades.  The  intersection  of  the  spectrum  with  u>  =  0  is  denoted  by  solid  black 
lines  to  better  illuminate  the  shape  of  the  spectrum 


From  Figs.  26b  and  26c  only  a  very  slight  asymmetry  is  noted  in  kz  which  arises  from  poorer  resolution 
in  this  direction  as  well  as  slight  tilting  of  the  PIV  cameras,  although  some  corrections  have  been  made 
for  the  latter.  For  all  practical  purposes,  these  spectra  are  symmetric  in  kz.  For  $(A:z,w)  this  symmetry 
arises  since  there  is  no  mean  spanwise  flow  and  thus  no  spanwise  directional  preference.  For  Q(kx,kz)  there 
should  physically  be  no  distinction  between  the  positive  and  negative  wavenumber  pairs.  Thus,  both  of  these 
spectra  can  be  represented  in  one  quadrant  with  the  magnitude  doubled  to  conserve  energy. 

For  Q(kx,uj)  the  two  quadrants  are  not  equivalent,  where  quadrant  I  is  interpreted  as  downstream  traveling 
waves  and  quadrant  II  as  upstream  traveling  waves.  The  amount  of  energy  in  quadrant  II  is  not  negligible 
(at  y+  =  108  quadrant  II  represents  about  12%  of  the  total  area  under  the  spectrum)  as  originally  suggested 
by  Morrison  and  Kronauer  [1969],  and  it  is  necessary  to  use  both  planes  when  integrating  Q(kx,u>)  to  recover 
cf>(kx),  or  any  other  integral  quantities,  otherwise  a  low  wavenumer  spectral  distortion  will  result. 

For  the  3D  spectrum  shown  in  Fig.  27,  the  same  symmetry  exists  in  kz  and  the  same  asymmetry  exists  in  kx 
as  might  be  expected  from  the  2D  spectral  plots.  Thus,  the  spectrum  can  be  represented  by  the  two  octants 
covering  all  kx  and  positive  kz  and  oj  with  the  magnitude  doubled  to  preserve  the  normalization  as  given  in 
Equation  27. 


JJJ  Q(kx,kz,uj)dkxdkz<Lj  =  A  JJJ  [©(fc^,  kz,u>)  +  Q(-kx,  kz,  w)]  dkxdkzduj 

—  oo  0 

oo 

4  JJj~  kxkzuj[@(kx,kz,uj)  +  @(-kx,kz,uj)]d\nkx  d\nkz  dlnuj 


(27) 


While  the  energy  in  the  upstream  traveling  waves  is  non-negligible  and  must  be  included  for  proper  normal¬ 
ization  of  Q(kx,  kz,uj ) ,  it  is  customary  in  the  literature  to  only  plot  data  for  positive  streamwise  wavenumber, 
kx,  when  these  spectra  are  presented  in  pre- multiplied  form.  While  this  may  seem  wrong  at  first,  when  pre¬ 
multiplied  and  presented  in  log  coordinates,  the  energy  in  the  first  decade,  which  encompasses  the  upstream 
traveling  waves,  is  spread  out  while  the  remaining  80%  of  the  energy  in  the  second  decade  is  concentrated 
and,  when  considering  absolute  energy  levels,  appears  much  more  energetic.  Therefore,  for  the  energy  lev¬ 
els  shown  in  the  following  figures,  the  upstream  traveling  waves  do  not  appear  so  plots  are  presented  over 
positive  kx  only. 


5.3.3  ID  Spectra 

Figs.  28a,  28c,  and  28e  show  the  ID  premultiplied  streamwise  velocity  spectrum  over  streamwise  wavenumber, 
spanwise  wavenumber,  and  frequency,  respectively,  for  the  uncorrected  data.  The  corrected  data  is  shown 
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(e)  u)<j>( u)  (f) 

Figure  28:  ID  premultiplied  streamwise  velocity  spectra  over  (a,b):  streamwise  wavenumber,  (c,d):  span- 
wise  wavenumbers,  and  (e,f):  frequency  where  the  lines  are  in  all  cases  — y+  =  34,  -  +  y+  =  108,  -  O 
y+  =  278.  The  uncorrected  data  is  shown  in  the  left  column  and  the  corrected  data  is  in  the  right  column. 
For  the  uncorrected  data,  the  dashed  lines  indicate  data  beyond  the  imposed  cutoff,  kcut ,  beyond  which  the 
data  begins  to  be  attenuated 
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Figure  29:  The  streamwise  wavenumber  and  frequency  spectra  are  compared  where  the  colors  are  black: 
kxfiikx),  gray:  Taylor  conversion  of  u><j){ui).  The  uncorrected  data  is  shown  in  the  left  and  the  corrected 
data  is  on  the  right.  For  the  uncorrected  data,  the  dashed  lines  indicate  data  beyond  the  imposed  cutoff, 
kCut  beyond  which  the  data  begins  to  be  attenuated.  Data  from  Ernr  and  Joubert  [1991]  at  Reg  =  1020 
is  represented  by  symbols  where  the  wall  normal  locations  are  (a,b):  y/5  =  0.04(+), 0.1(o),  (c,d):  y/5  = 
0.2(+), 0.35(0),  and  (e,f):  y/5  =  0.55(+) 
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Figure  30:  kxu>${kx,u>)  is  presented  for  all  three  planes.  The  uncorrected  data  is  presented  in  the  left 
column  and  the  corrected  data  is  in  the  right  column.  For  each  plot,  the  contours  represent  20%  to  80%  of 
the  maximum  energy  of  the  spectrum  moving  from  lighter  to  darker  shades  in  10%  increments 
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(a) 


(b) 


Figure  31:  (a):  The  convection  velocity  calculated  by  finding  the  line  of  local  maxima  along  kxu)&(kx,uj) 
at  y+  =  34  is  indicated  by  the  solid  black  line  (dashed  black  line  beyond  the  point  where  noise  influences  the 
calculation)  and  displayed  on  top  of  the  spectrum  where  the  contours  represent  20%  to  80%  of  the  maximum 
energy  of  the  spectrum  moving  from  light  to  dark  shades  in  10%  increments.  The  dash-dot  line  indicates 
where  u>  =  Ukx  signifying  a  convection  velocity  equal  to  the  local  mean.  The  four  open  circles  are  data 
from  Krogstad  et  al.  [1998]  at  the  same  y+  and  similar  Reynolds  number.  This  data  was  converted  from  its 
original  form,  convection  velocity  plotted  against  streamwise  probe  separation,  by  converting  the  separation 
to  a  wavenumber  and  using  the  convection  velocity  for  each  wavenumber  to  define  an  associated  frequency, 
(b):  The  convection  velocity  is  plotted  as  a  function  of  wavenumber  where  the  lines  and  symbols  are  the 
same  as  in  (a) 


in  Figs.  28b,  28d,  and  28f  for  comparison.  The  peak  in  the  spectrum  over  spanwise  wavenumber  moves 
to  lower  kz  (larger  spanwise  scales)  when  moving  away  from  the  wall  indicating  that  the  signature  of  the 
dominant  scales  is  growing  in  spanwise  extent  in  this  direction.  This  has  also  been  shown  in  the  channel 
flow  computations  of  Jimenez  et  al.  [2004]  and  boundary  layer  experiments  of  Tomkins  and  Adrian  [2005], 
among  others.  In  the  present  data,  as  the  peak  moves  to  higher  kz  closer  to  the  wall,  a  larger  energetic 
portion  of  the  spectrum  is  cut  off  by  spatial  averaging.  In  fact,  at  y+  =  34  a  better  resolution  would  be 
required  to  resolve  the  peak  in  the  spectrum  over  spanwise  wavenumbers.  Therefore,  any  integral  quantities 
in  the  z  direction,  including  the  rms  velocity  fluctuations  and  spectra,  will  be  more  in  error  as  the  wall  is 
approached.  This  is  likely  the  reason  for  the  slight  underestimation  of  urms  at  y+  =  34  shown  in  Fig.  25b. 

In  Fig.  29,  the  ID  spectra  over  frequency  and  over  streamwise  wavenumber  are  compared  to  one  another 
using  Taylor’s  hypothesis  and  to  the  hotwire  measurements  of  Erm  and  Joubert  [1991]  at  a  similar  Reynolds 
number  and  comparable  wall  normal  locations.  The  ID  spectra  are  in  fair  agreement  with  the  data  of  Erm 
and  Joubert  [1991],  and  this  agreement  is  improved  when  the  data  is  corrected  as  shown  in  the  plots  in  the 
right  column  of  Fig.  29.  This  improved  agreement  with  the  corrected  ID  spectra  should  also  be  indicative  of 
a  proper  correction  for  the  2D  and  3D  spectra  as  the  correction  is  applied  to  the  3D  spectrum  and  all  other 
spectra  are  calculated  via  integration  of  this  spectrum.  The  only  major  deviation  between  the  temporal  and 
spatial  spectra  is  at  y+  =  34.  This  deviation  may  be  indicative  of  the  variation  of  the  convection  velocity 
of  scales  from  the  local  mean,  although  noise  may  also  be  contributing  to  this  variation,  particularly  at  the 
high  kx  end,  as  shown  in  the  2D  spectra  in  Sect.  5.3.4.  Some  error  may  be  expected  at  this  wall  normal 
location  since  there  is  a  slight  deviation  of  the  mean  velocity  profile  from  the  measurements  of  DeGraaff  and 
Eaton  [2000]  as  shown  in  Figure  25a,  and  thus  other  statistics  may  be  adversely  affected.  For  the  other  two 
planes,  there  is  almost  no  deviation  between  the  spatial  and  Taylor  converted  temporal  spectra,  indicating 
that  the  frozen  flow  assumption  should  hold  at  these  locations. 

5.3.4  2D  Spectra 

The  2D  premultiplied  <&(kx,uj)  streamwise  velocity  spectra,  i.e.  kxoj&(kx,oj),  are  shown  for  each  plane  in 
Fig.  30  with  the  uncorrected  data  in  the  left  column  and  the  corrected  data  in  the  right  column.  Note  that 
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(e)  y+  =  278  (f)  y+  =  278 

Figure  32:  kzu><b(kz,uj)  is  presented  for  all  three  planes  where  the  contours  represent  20%  to  80%  of  the 
maximum  energy  of  the  spectrum  moving  from  light  to  dark  shades  in  10%  increments.  The  uncorrected 
data  is  presented  in  the  left  column  and  the  corrected  data  is  in  the  right  column.  For  the  uncorrected  data, 
the  unshaded  contours  are  regions  of  the  spectrum  beyond  the  streamwise  and  spanwise  cutoff,  kcut,  that 
may  be  irregularly  shaped  due  to  spectral  attenuation,  but  are  included  for  visualization 


53 


(e)  y+  =  278  (f)  y+  =  278 

Figure  33:  kxkz$(kx,kz)  is  presented  for  all  three  planes  where  the  contours  represent  20%  to  80%  of  the 
maximum  energy  of  the  spectrum  moving  from  light  to  dark  shades  in  10%  increments.  The  uncorrected 
data  is  presented  in  the  left  column  and  the  corrected  data  is  in  the  right  column.  For  the  uncorrected  data, 
the  unshaded  contours  are  regions  of  the  spectrum  beyond  the  streamwise  and  spanwise  cutoff,  kcut ,  that 
may  be  irregularly  shaped  due  to  spectral  attenuation,  but  are  included  for  visualization.  Linearly  spaced 
spectral  contours  from  del  Alamo  and  Jimenez  [2001]  for  a  channel  flow  computation  at  ReT  =  550  at  y+  =  90 
are  shown  with  colored  circles  where  cool  colors  indicate  less  energetic  portions  of  the  spectrum.  This  is 
compared  to  the  current  data  in  (c)  and  (d)  to  show  the  validity  of  the  correction 


54 


Figure  34:  kxkzu)Q(kx,kz,uj)  is  presented  at  each  wall  normal  location.  The  surfaces  represent  25%,  50%, 
and  75%  of  the  total  energy  in  each  spectrum  moving  from  light  to  dark  shades.  The  left  and  right  columns 
offer  two  different  views.  Data  beyond  the  kx  and  kz  cutoffs  are  not  shown 
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Figure  35:  The  contour  lines  represent  the  convection  velocity  at  y+  -  34  using  the  method  outlined  in 
Equation  28  where  solid  contours  denote  a  convection  velocity  equal  to  the  local  mean  and  dotted  contours 
with  symbols  are  convection  velocities  below  the  local  mean  with  o:  U  - uT,x :  U  -2ur,  and  ():  U  -3 uT. 
The  shaded  contours  represent  20  to  80%  of  the  maximum  energy  of  kxkzQ(kx,  kz )  in  10%  increments  moving 
from  light  to  the  dark  shades 


in  the  ID  spectra  of  Fig.  28,  the  maximum  energy  of  each  spectrum  decreases  as  measurements  are  taken 
further  from  the  wall,  as  expected.  In  the  2D  plots  here  and  3D  plots  to  follow,  contour  levels  represent 
a  fraction  of  the  maximum  energy  of  each  individual  spectrum  so  this  variation  will  not  be  present.  Also, 
note  that  the  total  energy  of  a  spectrum  is  different  before  and  after  correction,  so  the  levels  shown  for  the 
uncorrected  spectra  do  not  represent  the  same  levels  shown  for  the  corrected  spectra. 

One  notable  difference  between  the  corrected  and  uncorrected  kxu><&(kx,u>)  is  the  apparent  noise  that  appears 
in  the  corrected  spectrum  in  the  range  of  7  <  kx5  <  20  and  1  <  uiS/Uoo  <  9.  In  this  same  high  wavenumber 
range,  there  is  a  discrepancy  between  the  spatial  and  temporal  ID  spectrum  as  shown  in  Fig.  29a  and  more 
prevalently  in  the  corrected  spectrum  in  Fig.  29b  where  the  amplitude  of  the  spatial  spectrum  exceeds  the 
temporal  spectrum.  Since  the  ID  temporal  and  streamwise  wavenumber  spectra  should  agree  for  the  nearly 
homogeneous  small  scales,  the  discrepancy  beyond  kx5  =  7  likely  comes  from  noise.  It  is  hypothesized  that 
this  noise  in  the  near  wall  plane  is  from  the  measurement  of  nearly  stationary  particles  or  bubbles  that 
clustered  near  the  surface  and  were  illuminated  by  the  laser  sheet  during  PIV  acquisition. 

From  Q(kx,oj),  a  convection  velocity  can  be  calculated  for  each  streamwise  scale  as  shown  in  Fig.  31.  The 
method  of  finding  the  ridge  line,  the  line  of  maxima  along  the  spectrum  as  outlined  by  Goldschmidt  et  al. 
[1981],  is  used  here  where  the  convection  velocity  is  then  defined  as  ui/kx  along  the  ridge  line  (the  solid  line  in 
the  Fig.  31).  For  the  plane  shown  at  y+  -  34,  it  is  apparent  that  most  scales  travel  faster  than  the  local  mean 
velocity.  This  is  in  fair  agreement  with  Krogstad  et  al.  [1998],  whose  data  from  hotwire  measurements  is 
presented  for  comparison.  The  discrepancy  beyond  kxS  =  7  may  again  be  affected  by  the  noise  in  this  region 
as  described  previously,  and  for  this  reason  data  beyond  this  point  in  Fig.  31  has  been  denoted  by  a  dashed 
line.  In  addition,  the  slow  convection  velocities  for  low  wavenumbers  reflects  the  lack  of  scale  separation  at 
this  low  Reynolds  number,  namely  the  dominance  of  the  near-wall  structures  over  the  superstructures  which 
are  centered  at  locations  further  from  the  wall  but  inhabit  a  similar  wavenumber  range. 

Premultiplied  4 >(kz,u>)  and  <&(kx,kz)  are  shown  for  each  plane  in  Figs.  32  and  33,  respectively.  The  uncor¬ 
rected  data  is  shown  in  the  left  column  and  the  corrected  data  is  shown  in  the  right  column.  In  addition,  data 
from  the  channel  flow  computations  presented  by  del  Alamo  and  Jimenez  [2001]  at  y+  =  90  and  ReT  -  550  are 
included  for  comparison  to  the  current  data  at  y+  =  108  to  validate  the  correction  applied  (some  differences 
are  expected  due  to  the  difference  in  geometry  between  these  two  flows  as  well  as  the  slight  difference  in 
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Reynolds  number  and  wall  normal  location).  The  correction  appears  to  push  the  spectral  peak  in  the  right 
direction  and  align  the  current  data  to  the  data  from  del  Alamo  and  Jimenez  [2001].  The  dotted  contours 
for  the  uncorrected  data  show  the  region  that  is  beyond  the  wavenumber  cutoff  and  again  shows  that  the 
resolution  in  z  is  not  sufficient  to  resolve  the  peak  in  the  spanwise  direction  at  y+  =  34.  Again,  the  correction 
pushes  the  peak  in  the  spectra  to  higher  kz  and  increases  the  energy  in  the  spectra  in  general  near  this 
cutoff.  Upon  integration  in  kz,  the  resulting  corrected  ID  spectra  will  have  a  higher  level  throughout  than 
the  uncorrected  spectra,  as  evidenced  in  the  (j){kx )  and  4>(uj)  plots  in  Fig.  29. 

$>(kx,  kz )  shows  some  variability  between  the  three  planes  becoming  less  elongated  in  streamwise  wavenumber 
further  from  the  wall  as  shown  in  Figs.  33a,  33c,  and  33e.  The  range  of  energetic  streamwise  scales  narrows 
while  the  range  of  spanwise  scales  remains  fairly  constant  moving  further  from  the  wall  leading  to  a  more 
homogeneous  distribution  of  the  energy  in  the  wake  region. 

Note  that  all  2D  spectra  show  that  the  range  of  streamwise  scales  becomes  smaller  and  the  larger  streamwise 
scales  become  less  dominant  beyond  the  log  layer.  Recall,  that  these  conclusions  are  drawn  upon  low  Re 
data  so  the  scale  separation  is  not  large,  yet  this  conclusion  about  the  change  in  scale  size  from  the  inner 
layer  to  the  wake  region  should  apply  to  higher  Reynolds  number  data.  This  illustrates  the  difficulty  in 
performing  PIV  in  planes  near  the  wall  as  the  range  of  energetic  scales  is  broader  and  large  scales  play  a 
more  dominant  role.  Thus,  both  a  very  large  field  of  view  and  a  very  fine  resolution  are  needed,  which  was 
difficult  to  achieve  in  the  current  experiment.  The  best  solution  is  to  analyze  both  a  large  field  to  obtain 
the  large  scale  features  and  a  small,  well  resolved,  field  from  which  a  composite  spectrum  can  be  produced 
covering  the  whole  wavenumber  range,  which  is  the  focus  of  future  experiments. 

5.3.5  3D  Spectra 

Finally,  two  views  of  the  3D  premultiplied  streamwise  velocity  spectra,  kxkzu>@(kx,kz,uj),  are  presented  in 
each  wall  parallel  plane  in  Figs.  34a  to  34f  for  the  corrected  data.  The  uncorrected  data  has  been  omitted 
since  differences  between  the  two  are  already  clearly  outlined  by  the  2D  spectra  in  Sect.  5.3.4.  The  cutoff  of 
each  spectrum  is  clearly  shown  here,  and  the  effect  of  this  cutoff  on  all  of  the  spectra  presented  so  far  can 
easily  be  understood  by  considering  an  integration  in  one  or  more  coordinate  directions. 

These  three-dimensional  spectra  show  not  only  the  energy  distribution  over  all  scales,  but  over  all  scales 
traveling  at  all  velocities,  where  the  convection  velocity  of  a  scale  is  defined  as  uCtX  =  to/kx.  In  this  framework, 
if  the  flow  is  decomposed  into  traveling  waves  as  in  the  model  of  McKeon  and  Sharrna  [2010],  then  the 
spectra  are  the  footprints  of  these  waves  at  particular  wall  normal  locations.  For  flow  control  applications, 
measurements  taken  at  several  wall  normal  locations  would  provide  a  well  characterized  response  of  the 
boundary  layer  to  a  periodic  excitation  (as  studied  in  Jacobi  et  al.  [2010]  and  Jacobi  and  McKeon  [2011] 
using  hotwires)  and  alterations  of  the  spectra  could  be  used  to  see  how  a  particular  input  restructures  the 
flow.  Such  information  could  be  used  for  optimal  control  design. 

To  gain  insight  into  the  convection  velocity  of  different  scales  in  the  flow,  a  single  convection  velocity  can  be 
calculated  for  each  streamwise  and  spanwise  scale  pair  from  the  three-dimensional  spectra  using  equation  28 
from  del  Alamo  and  Jimenez  [2009]  and  used  in  earlier  works  by  Jimenez  et  al.  [2004]  and  Flores  and  Jimenez 
[2006]. 


uc{kx,kz,y ) 


1  f^0uQ(kx,kz,uj,y)duj 
kx  f^o®(kx,kz,u!,y)du> 


(28) 


This  convection  velocity  definition  uses  a  value  of  w  that  is  weighted  by  the  3D  streamwise  velocity  spectrum. 
In  this  way,  it  gives  the  dominant  convection  velocity  for  each  streamwise-spanwise  scale  pair.  A  map  of 
the  convection  velocity  is  presented  in  Fig.  35  where  lines  of  constant  convection  velocity  are  plotted  on 
top  of  kxkzQ(kx,kz)  to  highlight  the  convection  velocity  of  the  most  energetic  scales.  From  this  chart,  it 
is  apparent  that  most  scales  at  y+  =  34  travel  slower  than  the  local  mean  except  for  the  large  scales  in  the 
range  2  <  A x/5  <  5.  This  differs  from  the  trend  in  Fig.  31  where  almost  all  scales  recovered  traveled  faster 
than  the  local  mean.  This  difference  may  be  expected  as  the  method  of  calculating  the  convection  velocity 
differs  between  the  two  figures  where  one  searches  for  a  maximum  (method  for  Fig.  31)  while  the  other  looks 
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for  a  “center  of  mass”  (method  for  Fig.  35).  Also,  the  noise  that  exists  over  a  range  of  u>  as  well  as  the 
cutoff  imposed  in  u>  will  manifest  itself  differently  when  using  an  integral  method  such  as  Equation  28.  An 
investigation  of  such  convection  velocity  maps,  as  well  as  the  definition  of  convection  velocities  in  general,  is 
a  subject  of  future  work. 

5.4  Identification  of  swirling  coherent  structures 

5.4.1  SCS  identification  technique 

The  identification  of  swirling  coherent  structures  (SCS),  that  is,  coherent  structures  identified  by  regions  of 
swirling  flow  in  wall-parallel  planes  which  could  be  interpreted  as  wall-parallel  cuts  through  coherent,  three- 
dimensional  hairpin  vortices,  are  discussed  here.  Many  methods  have  been  proposed  for  identifying  vortices, 
and  subsequently  SCS,  but  as  discussed  by  Jeong  and  Hussain  [1995],  most  simple  intuitive  measures  fail  to 
identify  these  structures  correctly.  In  order  to  identify  vortex  cores  in  a  Galilean  invariant  setting,  methods 
including  the  A  [Chong  et  ah,  1990],  Q  [Hunt  et  ah,  1988],  A2  [Jeong  and  Hussain,  1995],  and  swirling 
strength  (A c,i)  [Zhou  et  al.,  1999]  criteria  have  been  proposed,  as  summarized  in  the  reviews  of  Jeong  and 
Hussain  [1995]  and  Chakraborty  et  al.  [2005].  These  methods  generally  provide  different  results  in  a  3D  flow 
field,  but  in  an  incompressible  flow  in  a  2D  plane,  as  studied  here,  all  of  the  criteria  give  identical  results. 

In  the  current  study,  the  wall-normal  component  of  the  swirling  strength,  given  in  Eqn.  29  for  a  2D2C  wall- 
parallel  measurement,  was  used  for  SCS  identification,  where  the  subscripts  x  and  z  represent  derivatives 
in  x  and  z,  respectively,  of  the  instantaneous  streamwise,  U,  and  spanwise,  W,  velocity.  We  follow  the 
convention  that  U  =  u  +  U,  (and  similarly  for  W)  where  U  is  the  instantaneous  velocity,  u  is  the  fluctuating 
velocity,  and  U  is  the  mean  velocity  averaged  in  the  homogeneous  spanwise  direction,  nearly  homogeneous 
streamwise  direction,  and  time,  and  thus  is  only  a  function  of  the  wall-normal  distance. 

Ac, i(x,  z )  =  | Im  [yJ{Ux  +  Wzy  +  4 (UXWZ  -  UZWX)]\  (29) 

A  signed  swirling  strength  was  also  used  to  identify  the  direction  of  swirl.  This  is  calculated  by  multiplying 
A c,j  by  the  sign  of  the  wall-normal  vorticity,  as  in  Natrajan  et  al.  [2007]. 

To  reduce  the  noise  in  the  numerical  derivatives  needed  to  calculate  A c,j,  the  velocity  field  was  convolved  with 
a  smoothing  kernel  before  differentiation.  Three  different  types  of  smoothing  were  investigated  including:  (1) 
A  moving  average  filter  with  no  weighting,  (2)  a  moving  average  filter  with  Gaussian  weighting,  and  (3)  the 
Savitzky-Golay  (SG)  filter  where  the  calculation  of  this  filter  and  the  smoothing  kernels  used  are  documented 
in  Savitzky  and  Golay  [1964].  To  compare  these  filters  and  chose  the  best  filter  size,  5  x  5,  7  x  7,  and  9x9 
smoothing  kernels  of  each  filter  type  were  tested  on  a  subsubsection  of  the  flow  field.  The  Gaussian  weighted 
average  as  well  as  the  SG  filter  performed  better  than  the  simple  moving  average  as  they  both  smoothed  the 
velocity  field  while  keeping  nearby  cores  separate.  Of  the  two,  the  SG  filter  (with  a  quadratic  polynomial) 
was  chosen  with  a  7  x  7  kernel  size,  as  it  allowed  derivatives  to  be  calculated  simply.  A  sample  velocity  field, 
with  the  signed  swirling  field  overlaid,  is  shown  in  Fig.  36.  Note  that  while  smoothing  is  necessary  to  reliably 
extract  SCS  from  the  velocity  fields,  the  accuracy  with  which  the  properties  of  the  SCS  can  be  measured  is 
reduced  as  the  signal  is  smeared  out  in  a  3Aa;  radius  (45 8V)  about  a  given  point  using  the  7x7  smoothing 
kernel. 

From  the  continuous  swirl  field,  individual  SCS  were  identified  by  extracting  regions  where  Ac,j  >  A rms,  the 
rms  swirling  strength,  which  can  be  considered  a  measure  of  background  noise.  This  threshold  was  also 
considered  by  Wu  and  Christensen  [2006],  although  they  found  a  more  stringent  threshold  of  Ac,j  >  1.5Arms 
separated  vortex  cores  better  for  their  measurements  in  the  wall-normal  plane.  For  the  sake  of  tracking 
SCS  (to  be  described  in  Subsection  5.5),  SCS  bounded  by  a  box  of  size  less  than  2  vector  spacings  in  the 
streamwise  or  spanwise  direction  were  neglected.  These  small  SCS  were  generally  found  to  be  energetically 
insignificant  and  would  often  disappear  from  one  timestep  to  the  next,  which  indicates  that  they  likely  arose 
from  noise  in  the  swirl  calculation. 

One  interpretation  of  persistent  SCS  is  the  wall-parallel  signature  of  circular  vortex  tubes,  perhaps  tilted 
with  respect  to  the  downstream  direction  and  possibly  consistent  with  the  footprint  of  hairpin  vortices.  In 
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Figure  36:  Clockwise  (red)  and  counter-clockwise  (blue)  swirling  motions  are  presented  on  top  of  a  velocity 
field  with  the  mean  velocity  subtracted,  shown  by  the  black  vectors,  at  y/S  =  0.014  at  one  point  in  time. 
Flow  is  from  left  to  right. 


Circular 
vortex  tube 


Figure  37:  The  SCS,  a  signature  of  a  circular  vortex  tube  cutting  the  wall-parallel  plane  of  view,  is  shown 
with  the  major,  a,  and  minor,  6,  axes  of  the  elliptical  SCS  labeled  along  with  the  vortex  tube  inclination 
angle,  a,  and  radius,  r. 


the  latter  case,  the  SCS  would  be  expected  to  consist  of  elliptical,  not  circular,  signatures  in  wall-parallel 
planes  as  discussed  by  Tomkins  and  Adrian  [2003]  and  illustrated  in  Fig.  37.  Under  this  interpretation, 
SCS  are  considered  to  be  the  signature  of  structures  which  are  coherent  in  space,  time  and  the  wall-normal 
direction,  and  hence  travel  at  a  y-independent  convection  velocity. 

5.4.2  SCS  identification  results 

From  each  extracted  SCS,  the  area,  A,  eccentricity  (of  the  ellipse),  e,  swirling  strength  weighted  centroid,  and 
average  core  velocity,  Ucore,  in  the  streamwise  direction  were  determined  (note  that  a  convection  velocity 
in  the  wall- normal  direction  cannot  be  determined  from  these  wall-parallel  measurements).  From  these 
measurements  the  equivalent  radius  of  the  SCS,  req,  and  the  inclination  angle  with  respect  to  the  streamwise 
direction,  a,  were  calculated.  Note  that  a  distinction  cannot  be  made  between  upstream  and  downstream 
inclination  from  our  data  since  the  elliptic  footprint  sketched  in  figure  37  would  appear  identical  in  the  wall- 
parallel  plane.  For  each  plane,  a  compilation  of  SCS  properties  in  500  statistically  independent  snapshots 
were  compared  to  the  literature  to  verify  the  identification  method. 
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Figure  38:  The  number  of  SCS  found  per  6x5  area  for  each  wall-normal  location 


The  average  number  of  SCS  identified  in  a  given  6x6  region  of  the  flow  is  shown  in  Fig.  38  for  each  wall- 
normal  location.  As  with  the  measurement  of  spanwise  vortex  heads  by  Wu  and  Christensen  [2006],  the 
current  data  also  shows  a  decrease  in  the  number  of  SCS  with  wall-normal  distance,  particularly  beyond  the 
log  layer,  consistent  with  the  attached  eddy  hypothesis. 

The  eccentricity  of  SCS,  e  =  yjl  -  ( b/a )2,  where  a  and  b  are  the  major  and  minor  axes  of  an  ellipse,  re¬ 
spectively,  was  calculated  by  constructing  an  ellipse  with  the  same  second  moments  of  area  as  the  identified 
swirling  region  and  measuring  the  ratio  of  the  major  and  minor  axes  of  this  construction.  From  the  eccentric¬ 
ity,  the  inclination  angle,  a,  of  a  circular  vortex  tube  cutting  a  wall-parallel  plane  (see  Fig.  38)  was  inferred 
using  Eqn.  30.  a  is  defined  as  the  angle  the  vortex  core  makes  with  respect  to  the  wall  and  ranges  from  0  to 
90°.  0°  indicates  a  vortex  core  locally  parallel  to  the  wall  (not  observable  in  the  present  experiments)  and 
90circ  indicates  a  vortex  core  locally  perpendicular  to  the  wall. 

a  =  sin-1  (\/l  -  e2)  (30) 

The  distribution  of  a  binned  to  the  nearest  5°  is  shown  in  Fig.  39a.  Note  that  the  pdf  presented  in  Fig.  39a, 
as  well  as  all  subsequent  pdfs  (unless  noted  otherwise),  was  normalized  such  that  the  area  under  the  pdf  was 
equal  to  1.  The  peak  near  40°  convincingly  demonstrates  that  the  extracted  SCS  are  elliptical  in  shape  on 
average.  In  addition,  the  a  distribution  does  not  vary  significantly  with  wall-normal  distance.  Comparison 
is  made  with  the  data  of  Ganapathisubramani  et  al.  [2006],  which  resolved  inclination  angles  ranging  from 
0  to  180°,  by  summing  the  portion  of  their  pdf  from  0  to  90°  with  that  from  180  to  90°  to  produce  a  pdf 
of  a  as  defined  in  the  present  work.  Further  investigation  by  Gao  et  al.  [2011],  however,  found  a  variation 
of  inclination  angle  (what  they  call  Qxz)  with  height  changing  from  25°  at  a  distance  of  y+  -  47,  to  45° 
at  y+  =  110,  to  50°  at  y+  -  198  for  a  similar  Reynolds  number.  Note  that  while  the  error  in  the  inferred 
inclination  angle,  which  arises  from  uncertainties  in  determining  the  length  and  width  of  the  SCS,  changes 
with  the  ratio  b/a  as  shown  in  Fig.  39b,  it  changes  little  when  b/a  <  0.5,  decreases  with  increasing  a,  and  is 
usually  below  about  20°  over  a  wide  range  of  a  and  b.  For  a  >  3 Ax,  the  error  falls  to  about  8°.  The  differences 
between  our  data  and  those  reported  by  Gao  et  al.  [2011]  close  to  the  wall  may  be  partially  explained  by 
the  differing  structure  detection  criteria,  and  perhaps  the  resolution  of  swirl  and  the  difficulty  identifying 
strongly  inclined  (low  a)  structures. 

Differences  between  the  present  data  and  the  data  of  Ganapathisubramani  et  al.  [2006]  in  Fig.  39a  can 
be  explained  as  follows.  First,  an  inclination  angle  of  0°  (indicating  swirling  motion  parallel  to  the  wall) 
cannot  be  inferred  from  the  present  2D2C  PIV  measurements.  Also,  for  all  core  sizes  considered,  there  were 
generally  2  or  fewer  inclination  angles  that  could  be  inferred  between  a  =  50°  and  a  =  90°  on  account  of  the 
the  discrete  pixel  locations,  with  the  number  of  angles  decreasing  with  decreasing  SCS  size,  as  indicated  by 
the  black  dots  in  Fig.  39b.  Thus  the  inferred  a  has  a  tendency  to  be  binned  to  50°  or  90°  degrees,  but  not  in 
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(b)  Error  in  angle  with  aspect  ratio 


Figure  39:  (a):  Distribution  of  the  SCS  inclination  angle,  a,  for  plane  1:  O,  2:  +,  3:  A,  and  4:  □.  Dashed 
black  line:  elevation  angle  distribution  from  Ganapathisubramani  et  al.  [2006]  (modified  as  described  in  the 
text)  at  y/S  =  0.09.  (b):  Error  in  the  inferred  inclination  angle,  Aa,  with  varying  aspect  ratio  for  an  elliptical 
SCS  with  a  major  axis,  a.  The  length,  a,  changes  from  3  to  6  vector  spacings,  increasing  in  the  direction  of 
the  arrow. 


between,  particularly  for  small  SCS.  The  deficit  of  counts  at  angles  beyond  50°  as  well  as  the  spike  at  90°, 
compared  to  the  pdf  from  Ganapathisubramani  et  al.  [2006],  likely  arose  from  this  poor  resolution  in  a. 

Since  the  SCS  were  shown  to  be  elliptical  in  wall-parallel  planes,  the  equivalent  radius,  req  -  \J  Ajn,  was  used 
as  a  measure  of  the  SCS  radius  to  validate  the  identification  method.  The  PIV  spatial  resolution  was  15 S„, 
hence  the  equivalent  diameter  would  only  be  measurable  to  this  value,  and  half  this  value  for  the  equivalent 
radius.  For  this  reason,  the  distribution  of  the  equivalent  radii  of  SCS  was  binned  to  the  nearest  7.2 8V 
as  shown  in  Fig.  40a.  As  noted  previously,  though,  filtering  the  velocity  field  before  calculating  the  swirl 
field  decreased  the  accuracy  with  which  the  vortex  radius  could  be  measured,  with  the  smallest  measurable 
radius  15 5V,  potentially  smeared  out  to  60<5„.  The  change  in  the  mean  radius  with  wall-normal  location  is 
recorded  in  Fig.  40b.  While  the  two  planes  near  the  wall  share  an  almost  identical  radius  distribution,  as 
the  free  stream  is  approached,  the  distribution  begins  to  skew  towards  larger  core  radii  consistent  with  the 
wall-normal  growth  in  eddy  dimension  associated  with  Townsend’s  attached  eddy  hypothesis.  The  lack  or 
the  variation  of  the  core  radii  for  planes  nearer  the  wall  may  be  at  least  partially  explained  as  an  artifact  of 
the  smoothing  of  the  swirl  field  already  discussed. 

Comparing  to  the  literature,  Carlier  and  Stanislas  [2005]  measured  a  mean  radius  of  20-30<5„,  increasing 
with  wall-normal  distance.  Differences  between  the  current  measurements  and  Carlier  and  Stanislas  [2005] 
can  be  rectified  by  calculating  the  radius  of  the  vortex  tube  cutting  the  wall-parallel  plane,  inferred  from 
the  inclination  angle  and  the  aspect  ratio  of  the  SCS  measured  in  plane.  Making  reference  to  Fig.  37  and 
using  the  relationship  req  =  \J Ajn,  the  area  of  an  ellipse  A  =  nab,  assuming  a  tilted  circular  vortex  tube 
a  =  b/sin(a),  and  the  core  radius,  r  =  b,  then: 


r  =  req\J  sin(a)  (31) 

By  applying  this  correction,  the  corrected  radii  fall  within  the  expected  range  of  sizes  as  shown  by  the  triangle 
symbols  in  Fig.  40b.  Thus  our  identification  method  gave  results  consistent  with  the  existing  literature  based 
on  ensemble  averages  of  spatial  fields  and  an  analysis  of  the  temporal  behavior  of  the  SCS  was  performed. 
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(a)  Core  radius  distribution  (b)  Mean  core  radius 

Figure  40:  (a):  Distribution  of  the  SCS  equivalent  radii,  req ,  for  plane  1:  O,  2:  +,  3:  A,  4:  □.  (b): 
Wall-normal  variation  of  the  mean  value  of  O:  reqi  A:  r  (calculated  from  req  using  Eqn.  31),  and  x:  the 
vortex  radius  from  the  wall-normal  measurements  of  Carlier  and  Stanislas  [2005]  at  Reg  =  7500. 


5.5  Swirling  coherent  structure  tracking 

5.5.1  SCS  tracking  technique 

The  swirl-weighted  centroids  of  SCS  were  tracked  in  time  using  the  time-resolved  data  to  analyze  SCS 
trajectories  in  the  wall-parallel  plane  as  well  as  to  investigate  their  average  lifetime  and  properties.  Tracking 
was  performed  on  Eulerian  structures;  that  is  structures  were  identified  in  each  frame  independent  of  one 
another  as  opposed  to  tracking  the  evolution  of  Lagrangian  tracers  in  time,  e.g.  Yang  and  Pullin  [2011].  The 
velocity  fields  were  subsampled  in  time  so  that,  on  average,  a  SCS  would  travel  forward  at  least  one  vector 
spacing  between  frames;  A tnew  =  A x/U(y). 

The  SCS  centroids  in  consecutive  frames  were  compared  to  determine  the  SCS  (centroid)  trajectories.  For 
each  SCS  in  the  first  frame,  a  nearest  neighbor  with  the  same  swirling  direction  was  identified  in  the  second 
frame.  If  the  nearest  neighbor  was  within  a  “reasonable”  distance  of  the  original  SCS,  then  the  two  were 
identified  as  being  the  same  structure.  The  criteria  for  a  valid  nearest  neighbor  required  that  the  SCS  either 
not  move  in  the  streamwise  direction  (for  any  convection  slower  than  the  local  mean,  the  SCS  will  not  move 
a  whole  vector  spacing)  or  that  it  did  not  move  faster  than  the  free-stream  velocity.  This  range,  defined  as 
0  <  (aq  -  Xq)  <  Uoo/U(y)Ax,  was  widened  to  account  for  inaccuracies  in  determining  the  SCS  centroid  by 
Aa;  on  both  ends.  In  addition,  spanwise  motion  was  restricted  such  that  \z2  -  z i|  <  2 A z  per  time  step.  This 
constraint  was  applied  to  assure  that  neighboring  SCS  were  not  mistakenly  picked  up  in  the  tracking  method 
and  was  deemed  a  reasonable  limit  given  that  the  main  direction  of  flow  was  in  the  streamwise  direction  and 
the  spanwise  motion  was  observed  to  be  only  a  small  fraction  of  this. 

Regardless  of  the  restrictions  placed,  there  were  still  issues  with  tracking  SCS  in  time;  one  example  is  shown 
in  Fig.  41.  In  the  first  frame,  two  nearby  SCS  of  the  same  sign  are  illustrated  by  blue  circles  with  the  x 
marking  their  centroids.  In  the  second  frame,  the  SCS  are  wrongly  marked  as  a  single  core  due  to  poor 
spatial  resolution  of  the  nearby  structures,  and  thus  only  one  centroid  is  found.  The  circled  SCS  in  the  first 
frame  has  now  lost  its  trace  as  only  one  of  the  two  SCS  can  be  linked  with  this  new  structure.  In  the  final 
frame,  the  SCS  are  again  correctly  identified,  but  the  circled  SCS  is  incorrectly  marked  as  a  new  SCS  as 
there  was  no  distinct  nearest  neighbor  in  the  previous  frame.  To  help  alleviate  this  issue,  after  identifying 
all  SCS  trajectories,  the  SCS  position  at  the  end  of  a  trajectory  was  projected  forward  at  the  local  mean 
velocity  for  the  next  two  frames.  Following  Fig.  41,  the  position  of  the  circled  SCS  in  frame  1  would  be 
projected  forward  in  frames  2  and  3.  Then,  if  the  start  of  a  new  SCS  trajectory  was  found  to  be  a  nearest 
neighbor  of  that  projected  SCS,  the  two  trajectories  would  be  merged  together.  Referring  to  Fig.  41,  the 
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Figure  41:  Cartoon  showing  how  a  SCS  trace  can  be  lost  due  to  poor  resolution.  SCS  are  represented  by 
blue  regions  and  the  SCS  centroids  are  marked  with  red  x’s.  The  dashed  circles  indicate  where  tracking  of 
a  particular  SCS  has  failed. 
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Figure  42:  SCS  trajectories  at  y/S  =  0.3  that  persist  for  more  than  one  eddy  turnover  time  (S/Uoo) 


SCS  trajectory  ending  with  the  circled  SCS  in  frame  1  would  be  merged  with  the  trajectory  staring  with  the 
circled  SCS  in  frame  3. 


5.5.2  SCS  tracking  results 

The  SCS  were  found  to  move  predominantly  in  the  streamwise  direction,  as  expected,  with  limited  motion 
in  the  spanwise  direction,  as  shown  in  the  trajectories  of  Fig.  42.  Note  that  no  additional  smoothing  has 
been  applied  to  these  trajectories  in  order  to  maintain  a  faithful  representation  of  the  spanwise  resolution. 
The  average  trajectory  angle,  0t,  was  calculated  by  fitting  a  straight  line  to  each  trajectory.  A  histogram  of 
the  trajectory  angles  is  shown  in  Fig.  43a.  The  width  of  these  histograms,  and  thus  any  spanwise  motion,  is 
small,  and  decreases  with  wall-normal  distance.  At  all  wall-normal  locations,  one  would  expect  a  trajectory 
angle  of  0°  on  average  due  to  the  symmetry  of  the  pdf  about  0°. 

The  dominant  peak  in  the  6t  histograms  at  0°,  along  with  observations  from  the  time-resolved  velocity 
fields,  suggests  that  vortex  packets  on  average  travel  on  nearly  straight  trajectories,  deviating  very  little  in 
their  spanwise  location  over  the  entire  field  of  view,  but  likely  growing  in  the  vertical  direction.  It  seems 
unlikely  that  the  well-known  spanwise  growth  with  distance  from  the  wall  occurs  by  continuous  expansion  of 
individual  SCS  in  the  spanwise  direction,  in  time,  as  they  convect  downstream,  but  rather  through  a  discrete 
mechanism  such  as  the  packet  merging  mechanism  of  individual  component  vortices  discussed  by  Tomkins 
and  Adrian  [2003].  Note,  however,  that  the  combination  of  our  tracking  algorithm  and  spatial  resolution  is 
not  sufficiently  sophisticated  to  distinguish  convincingly  between  such  mechanisms. 

Using  the  time-resolved  data,  the  distribution  of  the  SCS  lifespans  was  determined,  where  TL  denotes  the 
lifespan  of  a  SCS.  The  lifespan  is  defined  as  the  amount  of  time  that  a  SCS  in  one  of  the  wall-parallel  planes 
can  be  tracked  in  time.  The  lifespan  is  skewed  toward  lower  values  due  to  SCS  which  leave  the  field  of  view  in 
x  or  z.  Also,  it  is  possible  for  a  SCS  to  travel  up  or  down  in  the  wall-normal  direction  (y)  as  time  progresses 
or  to  grow  away  from  the  wall  with  the  boundary  layer,  which  grows  by  about  10%  over  the  field  of  view, 
and  thus  leave  the  field  of  view  and  further  decrease  the  measured  lifespan.  With  these  constraints  in  mind, 
a  pdf  of  the  SCS  lifespan  for  each  plane  is  shown  in  Fig.  43b.  A  majority  of  the  SCS  persist  for  less  than 
1  eddy  turnover  time,  Te,  with  an  initial  decay  rate  of  1.6  decades  per  Te  for  all  wall-normal  locations.  In 
all  cases,  the  pdf  has  dropped  to  0.01-0.02  by  Te  -  Te-  It  must  be  noted  that  most  of  these  counts  at  small 
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(a)  Trajectory  angle 


Figure  43:  (a):  Histogram  of  the  SCS  trajectory  angles  for  planes  1:  O,  2:  +,  3:  A  and  4:  □.  Data  is 
compiled  from  both  CW  and  CCW  swirling  motions  and  for  all  runs  recorded,  (b):  Pdf  of  SCS  lifespan  for 
planes  1:  O,  2:  +,  3:  A,  and  4:  □.  Dashed  line:  the  slope  of  all  pdfs  before  0.56/Uoo 


Tl  are  likely  due  to  issues  with  tracking  a  particular  SCS  in  time,  such  as  the  issue  illustrated  in  in  Fig.  41. 
As  the  free  stream  is  approached,  SCS  tend  to  live  for  longer  periods  of  time,  as  noted  by  the  larger  values 
of  the  pdf  at  long  times  for  planes  3  and  4  in  Fig.  43b,  which  is  postulated  to  be  an  effect  of  the  reduction 
of  turbulent  stresses  deforming  and  possibly  destroying  these  structures  further  from  the  wall.  The  plateau 
in  the  pdf  at  long  times  (between  3.5  and  5  TeUoo/5)  in  plane  4  is  an  effect  of  the  limited  observation  time 
(limited  by  the  size  of  the  domain)  as  will  be  discussed  later  in  this  subsection  in  conjunction  with  Fig.  46. 

A  simple  estimate  of  the  SCS  convection  velocity  was  extracted  from  the  SCS  trajectories  by  dividing 
the  change  in  streamwise  position  of  the  SCS  centroid  by  the  time  between  the  two  frames  investigated. 
Note  that  this  linear  scheme  cannot  identify  changes  in  streamwise  convection  velocity,  which  have  been 
observed  to  occur  over  timescales  larger  than  the  Lagrangian  timescale  by  Elsinga  et  al.  [2012],  spanwise 
convection  velocity,  or  a  non-zero  wall-normal  convection  velocity,  but  does  provide  a  robust  characterization 
of  convection  velocities  over  our  large  wall-parallel  domain.  Of  course  more  sophisticated  schemes  could  be 
applied.  Between  individual  time  steps,  the  measured  convection  velocity,  Utraj  fluctuated  due  to  issues  in 
finding  the  SCS  centroid.  In  particular,  expecting  a  positioning  error  of  half  the  window  length,  7.56 v,  would 
result  in  an  error  in  the  measured  velocity  of  u*rror  =  7.1  -9.6.  Over  a  number  of  time  steps  the  average 
convection  velocity  quickly  converged  to  a  mean  value  and  this  error  shrinks.  For  this  reason,  the  average 
convection  velocities  over  a  trajectory,  <  Utraj  >  and  <  Utraj  >  -  <  Utraj  >  —U  are  reported,  where  the  triangle 
brackets  indicate  an  average  over  the  entire  vortex  trajectory. 

The  pdf  of  <  Utraj  >  for  SCS  whose  lifespan  exceeded  one  eddy  turnover  time,  i.e. ,  Tf,  >  Te,  is  shown  in 
Fig.  44a  for  each  wall-normal  location.  For  these  cases,  the  convection  velocity  has  been  calculated  and 
averaged  over  at  least  16  time  steps  ( t+  >  5.44).  The  peaks  of  all  the  pdfs  are  close  to  the  local  mean,  U(y), 
with  the  peak  shifting  toward  slightly  slower  values  as  y  increases.  All  pdfs  also  have  a  similar  low-speed  tail 
which  is  longer  than  the  high-speed  tail.  The  major  differences  between  different  wall-normal  locations  is  a 
thinning  of  the  high-velocity  end  of  the  pdf,  and  hence  an  increase  in  the  peak  amplitude  of  the  pdf,  most 
notable  at  y/5  =  0.48,  an  effect  of  the  upper  limit  on  the  convection  velocity  set  by  the  free-stream  velocity. 
This  is  readily  observable  in  the  pdf  of  <  Utraj  >  in  Fig.  44b  and  is  also  evident  in  the  estimates  of  the 
convection  velocity  of  spanwise  vortex  cores  obtained  by  integration  of  the  spatial  velocity  field  by  Wu  and 
Christensen  [2006].  It  is  apparent  that  no  SCS  traveled  at  speeds  above  the  free-stream  velocity  (U£,  -  22.5) 
and  none  were  observed  below  U+  =  8,  the  estimate  of  the  lower  bound  on  the  convection  velocity  of  energetic 
disturbances  in  the  boundary  layer  identified  by  Morrison  et  al.  [1971]. 

The  average  of  all  velocity  vectors  within  each  SCS,  Ucore  (a  commonly-used  approximation  to  the  convection 
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Figure  44:  (a):  Pdf  of  <  Utraj  >  for  planes  1:  O,  2:  +,  3:  A,  and  4:  □.  (b):  Distribution  of  SCS  trajectory 
velocity,  <  Utraj  >  for  plane  1:  O,  2:  +,  3:  A,  and  4:  □. 


velocity  for  non-time-resolved  PIV  studies  -  note  that  Utraj  cannot  be  calculated  from  an  ensemble  of 
instantaneous  fields),  was  also  calculated  as  a  measure  of  the  SCS  convection  velocity.  It  is  equivalent  to 
the  convection  velocity  that,  when  subtracted  from  the  instantaneous  velocity,  would  produce  a  set  of  closed 
streamlines  at  the  location  of  the  SCS.  For  comparison,  the  calculated  values  are  averaged  over  the  vortex 
trajectory  to  obtain  <  Ucore  >  =  <  ucore  >  +U.  The  mean  value  of  <  ucore  >  over  all  trajectories,  <  ucore  >m, 
is  compared  to  <  Utraj  >m  at  each  wall- normal  location  in  Fig.  45a,  showing  values  of  <  ucore  >m  closer  to 
the  local  mean  velocity  than  <  utraj  >mi  with  a  value  slightly  higher  than  the  local  mean  near  the  wall  and 
slower  than  the  mean  far  from  the  wall.  The  finding  that  <  Utraj  >m  is  lower  than  the  local  mean  at  all 
wall-normal  distances  is  similar  to  the  findings  of  Gao  et  al.  [2011]  that  vortices  propagate  at  96-98%  of  the 
local  mean  with  stronger  vortices  convecting  slower  than  weaker  ones. 

To  explain  the  difference  between  <  itcore  >m  and  <  Utraj  >m,  the  convection  velocity  pdfs  of  each  at  plane 
3  are  compared  in  Fig.  45b.  This  figure  shows  a  similar  distribution  of  velocities  above  the  local  mean, 
but  the  <  Utraj  >  pdf  has  a  longer  low  speed  tail  and  lower  central  peak  than  the  <  ucore  >  pdf,  and  hence 
a  lower  mean  value.  While  the  distribution  of  convection  velocities  is  similar  for  <  ucore  >  and  <  Utraj  > 
at  all  wall-normal  locations,  the  Utraj  measure  identifies  more  slow-moving  SCS.  This  effect  could  be  at 
least  in  part  attributable  to  the  importance  of  the  (unquantified)  wall- normal  convection  velocity,  Vtraj ,  of 
the  SCS.  Under  the  assumption  of  a  <  vtraj  >  pdf  with  a  width  of  the  same  order  as  for  the  streamwise 
convection  velocity,  the  error  associated  with  the  neglected  wall-normal  motion  can  be  determined  to  lead  to 
a  broadening  of  the  <  utraj  >  pdf  that  is  not  negligible  with  respect  to  the  differences  identified  in  Fig.  lib. 
Positive  <  Vtraj  >  could  give  rise  to  longer  apparent  tails  of  the  <  Utraj  >  pdf.  The  temporal  evolution  and 
growth  of  SCS-like  structures  is  well-known.  However,  in  the  absence  of  synchronous  measurements  of  the 
PIV  measurements  in  the  wall-normal  direction,  no  conclusive  statement  can  be  made  about  the  magnitude 
and  importance  of  any  wall-normal  convection  velocity.  It  is  believed  that  the  trends  discussed  below  are 
robust  to  this  potential  source  of  error. 

Finally,  to  complete  the  picture,  the  lifespan  and  convection  velocity  of  SCS  are  considered  simultaneously 
using  joint  pdfs  as  presented  in  Fig.  46  for  each  wall-parallel  plane.  The  red  line  indicates  the  edge  of  the 
domain  determined  by  dividing  the  size  of  the  streamwise  field  of  view  by  the  convection  velocity,  <  Utraj  >■ 
This  edge  is  the  time  at  which  a  SCS  would  leave  the  field  of  view  if  the  SCS  started  at  the  most  upstream 
end  of  the  field  of  view  and  was  traveling  at  a  fixed  velocity,  <  >.  A  large  buildup  of  counts  at  this 

limit  is  found  at  y/S  =  0.48  as  shown  in  Fig.  46d.  Thus,  in  this  plane,  the  lifespan  of  a  reasonable  portion  of 
the  SCS  persist  for  longer  than  5  eddy  turnover  times.  This  buildup  of  counts  was  previously  identified  as 
a  plateau  in  the  SCS  lifespan  histograms  shown  in  Fig.  43b  (i.e.,  the  integral  of  the  jpdf  over  <  ■  >)  and 
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Figure  45:  (a):  Mean  convection  velocity  for  each  plane  where  □:  <  'utrQ)  >m  and  O:  <  ucore  >m.  (b): 
Comparison  of  convection  velocity  pdfs  for  plane  3.  Dashed  line:  <  ucore  >,  solid  line:  <  Utraj  > 


confirms  the  trend  of  increasing  lifespan  with  wall-normal  distance.  Despite  the  relatively  large  field  of  view 
in  this  study,  it  still  places  a  constraint  on  the  statistics  associated  with  long  life-time  SCS. 

5.6  Implications  for  models  of  swirling  coherent  structure 

The  subsections  above  present  the  statistical  features  associated  with  swirling  coherent  structures  identified 
non-synchronously  in  four  wall-parallel  planes.  In  what  follows,  we  exploit  converged  statistics  over  many 
realizations  of  SCS,  rather  than  a  rigorous  connection  of  instantaneous  structures,  to  investigate  how  well 
the  data  agree  with  a  model  of  SCS  behavior  that  constrains  the  structures  to  travel  at  or  close  to  an 
identical  convection  velocity  at  all  wall-normal  planes  in  order  to  remain  coherent  for  a  significant  time. 
We  then  ascribe  observations  in  different  planes  to  distinct  “model  structure  types”  with  similar  convection 
velocities.  Of  course  the  decorrelation  effect  associated  with  the  full  turbulent  velocity  field  may  prevent  a 
single  structure  from  being  observed  to  span  multiple  planes,  however  structures  traveling  with  the  same 
convection  velocity  will  be  understood  to  have  a  similar  origin. 

Such  a  model  is  entirely  consistent  with  the  attached  eddy  ideas  proposed  by  Townsend  [1976]  and  developed 
further  by  Perry  and  co-workers,  e.g.  Perry  and  Chong  [1982],  Perry  et  al.  [1986],  Perry  and  Marusic  [1995]. 
These  authors  developed  a  model  for  the  underlying  mechanics  of  wall  turbulence  based  on  the  action  of  self¬ 
similar  hierarchies  of  hairpin-like  vortices  that  scale  on  distance  from  the  wall.  The  wall-normal  extent  of  each 
such  vortex  would  be  understood  to  be  associated  with  a  single  structure  convection  velocity.  The  subsequent 
conceptual  model  of  vortical  structure  developed  by  Adrian  et  al.  [2000],  illustrated  in  their  Fig.  25,  outlines 
packets  of  hairpin  vortices  that  travel  with  similar  convection  velocities.  By  extension,  for  such  packets  to 
have  temporal  coherence,  individual  structural  components  must  be  assumed  to  convect  at  similar  speeds 
across  their  full  wall-normal  extents.  This  assumption  is  supported  by  the  film  of  Bandyopadhyay  and  Head 
[1979],  which  shows  that  all  points  on  a  given  vortex,  identified  using  smoke  visualization,  appear  to  convect 
at  a  uniform  streamwise  velocity.  However,  indications  from  DNS  [Lozano-Duran  et  al.,  2012]  suggest  that 
wall-normal  variation  of  the  convection  velocity  for  a  range  of  structures  is  likely,  leading  necessarily  to 
reduced  temporal  coherence  as  structures  are  stretched  and  eventually  torn  apart. 

In  what  follows,  we  investigate  model  structures  whose  signature  in  different  planes  is  constrained  to  travel 
at  the  same  convection  velocity.  While  multi-plane,  synchronous  PIV  measurements  would  be  required  to 
prove  the  hypothesized  relationships,  we  believe  that  our  approach  suggests  a  useful  and  statistically-sound 
means  of  classifying  possible  SCS  connectivity  for  the  cases  where  individual,  constant  structure  convection 
velocities  are  appropriate.  We  note  that  the  discussion  that  follows  is  limited  to  these  cases  only. 
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(c)  Plane  1,  y+  =  122 


(d)  Plane  1,  y+  =  198 


Figure  46:  The  joint  pdfs  of  SCS  convection  velocity  and  SCS  lifespan  for  plane  (a):  1,  (b):  2,  (c):  3,  and 
(d):  4.  The  color  indicates  the  probability  of  having  a  vortex  with  a  given  lifespan  and  trajectory  velocity, 
where  the  double  integral  over  <  U^ra]  >  and  T^Um/S  is  equal  to  1.  Black  line:  local  mean  velocity,  dashed 
red  line:  edge  of  the  domain 
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(a)  Histogram  comparison  (b)  Cartoon  of  inferred  structures 

Figure  47:  (a)  Comparison  of  convection  velocity  histograms  at  y/6  =  0.08  and  y/6  =  0.30  leading  to  the 
identification  of  3  distinct  regions.  The  labeled  regions  are  referenced  in  (b)  showing  the  structure  shapes 
inferred  from  the  histograms  in  (a). 


Fig.  47  shows  an  overlay  of  convection  velocity  histograms  for  two  different  planes.  Region  B,  where  the 
histograms  overlap,  represents  structures  that  convect  at  the  same  velocity  and  could  be  interpreted  to  be 
connected  between  the  two  planes.  Regions  A  and  C  represent  structures  that  exist  in  only  one  of  the  two 
planes.  By  analyzing  these  comparisons  for  all  combinations  of  the  four  measured  planes,  ten  different  model 
structures  types  can  be  inferred  as  shown  in  Fig.  48,  organized  by  likelihood  of  observation.  The  percentage 
of  each  type  of  structure  in  the  flow  is  also  recorded  in  this  figure  (identified  from  the  four  curves  in  Fig.  44b). 
To  label  the  structures,  we  introduce  the  terminology  where  i  and  j  are  the  numbers  of  the  plane 

nearest  the  wall  and  furthest  from  the  wall  through  which  the  structure  passes,  respectively.  The  method  for 
counting  and  calculating  the  convection  velocity  histograms  for  these  structures  is  described  in  Appendix  ??. 

Results  from  this  analysis  are  described  in  the  remainder  of  this  subsection,  but  first,  some  limitations  of 
this  method  are  outlined  or  reiterated.  The  first  issue  to  note  is  that  the  connectivity  between  wall-parallel 
planes  depends  on  the  threshold  chosen  to  identify  SCS.  For  example,  with  a  lower  threshold,  more  SCS 
would  be  identified  in  all  planes,  and  this  could  lead  to  an  increase  in  the  connectivity  between  different 
planes,  possibly  increasing  the  number  of  SCS  that  extend  to  planes  1  and  4.  It  has  also  been  assumed  that 
the  convection  velocity  is  y-independent  and  constant  in  time,  which  does  not  allow  for  the  growth  of  SCS, 
that  is  they  cannot  stretch  in  the  streamwise  direction  as  they  travel  downstream.  Similarly,  this  framework 
does  not  have  a  way  to  count  vortices  which  grow  in  wall-normal  extent  and  reach  a  plane  further  from  the 
wall  at  some  later  point  in  time. 

With  these  constraints  in  mind,  some  statements  can  be  made  about  the  likely  distributions  of  model 
structure  type.  From  Fig.  48  it  is  apparent  that  the  most  prevalent  SCS  is  ^(l^),  a  structure  that  reaches 
down  to  plane  1  ( y+  =  38)  and  does  not  reach  beyond  the  log  layer  (plane  2).  We  term  this  an  attached 
structure,  implying  only  that  the  footprint  of  the  structure  reaches  down  to  the  wall.  The  taller  attached 
structures,  ^(i^)  and  4'(1,4),  are  less  common  than  ^(l^)  (which  occurs  in  22%  of  all  observations), 
in  agreement  with  the  observation  that  the  majority  of  hairpin  heads  are  found  in  the  log  layer  [Wu  and 
Christensen,  2006].  Fig.  49a  shows  the  variation  of  the  number  of  attached  structures  with  wall-normal 
distance,  an  approximately  logarithmic  decay.  This  decrease  in  attached  structures  with  wall-normal  extent 
is  supported  by  the  measurements  Hutchins  et  al.  [2005],  where  streamwise  velocity  correlations  were  used 
to  infer  that  a  majority  of  structures  should  be  detached  by  y/S  =  0.4  and  at  least  some  should  be  detached 
as  close  as  y/S  =  0.05. 

Jpdfs,  as  shown  in  Fig.  46,  can  also  be  constructed  for  each  model  structure  type.  Similar  to  the  process 
described  in  Fig.  47,  joint  lifespan,  convection  velocity  histogram  surfaces  are  constructed  for  two  different 
planes  and  overlaid.  Regions  enclosed  by  both  histogram  surfaces  can  be  interpreted  as  model  structures 


68 


m4) 


Figure  48:  The  model  structures  are  shown  in  order  of  their  prevalence  in  the  boundary  layer.  The 
percentage  given  is  rounded  to  the  nearest  percent  and  SCS  with  similar  percentages  are  grouped  together. 
Lighter  colored  SCS  are  attached  and  darker  colored  SCS  are  detached  from  the  wall. 


Figure  49:  (a):  The  percentage  of  SCS  at  a  given  wall-normal  distance  that  are  “attached”  to  the  plane 
nearest  the  wall,  (b):  Lifespan  pdfs  for  all  attached  structures  (A)  and  all  detached  structures  (□).  Dashed 
black  lines  (no  symbols):  log  fit  to  the  lifespan  decay  for  attached  or  detached  structures  up  to  the  vertical 
dashed  line  with  the  corresponding  symbols 
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(a)  Plane  2 


(b)  Plane  4 


Figure  50:  The  trajectory  velocity  pdf  for  planes  2  and  4  are  shown  by  a  solid  line  in  (a)  and  (b),  respectively 
and  the  vertical  dashed  lines  indicate  the  local  mean  velocity  in  plane  1:  O,  2:  +,  3:  A,  and  4:  □.  The 
other  lines  are  the  convection  velocity  pdfs  for  (a):  attached  structures  passing  through  plane  2;  4/(1, 2):  +, 
4/(1, 3):  A  and  4/(1, 4):  □,  (b):  detached  structures  passing  through  plane  4;  4/(2, 4):  +,  4>(3,4):  A,  and 
4/(4, 4):  □. 


connected  between  the  two  planes,  and  regions  enclosed  by  only  one  of  the  two  histograms  as  structures 
that  exist  in  only  one  of  the  two  planes.  The  process  described  in  Appendix  ??  can  also  be  applied  to  these 
joint-histograms  to  determine  the  distribution  4 /(*,_))  for  each  structure  type. 

By  integrating  the  jpdf  for  a  particular  model  structure  type  over  the  trajectory  velocity,  the  lifespan  of 
that  structure  type  can  be  found.  The  lifespan  pdfs  for  all  attached  structures  were  similar  as  well  as  all 
detached  structures,  so  an  average  lifespan  pdf  is  presented  for  attached  and  detached  structures  in  Fig.  49b. 
The  faster  decay  rate  of  the  attached  structures  compared  to  the  detached  structures  is  notable,  with  the 
probability  of  having  a  SCS  exist  for  an  additional  eddy  turnover  time  dropping  by  less  than  half  a  decade 
(0.4)  for  detached  structures,  but  dropping  by  nearly  a  full  decade  for  attached  structures.  The  pdf  at  large 
Tl  is  less  converged  for  attached  SCS  due  to  the  small  number  of  counts  at  lifespans  beyond  4 Te  and  the 
decay  at  large  times  for  the  detached  structures  is  due  to  the  limited  field  of  view  beyond  about  5 Te  as 
noted  earlier. 

The  convection  velocity  pdfs  associated  with  each  type  of  model  structure  can  be  calculated  by  integrating 
the  jpdfs  over  the  lifespan.  The  velocity  pdfs  of  different  structure  types  can  be  used  to  elucidate  the  makeup 
of  the  convection  velocity  pdfs  presented  in  Fig.  44.  Two  examples  are  examined  in  Fig.  50.  The  convection 
velocity  pdfs  of  all  attached  and  all  detached  structure  types  are  shown  in  Figs.  50a  and  50b,  respectively. 
These  structure  pdfs  have  been  normalized  considering  the  percentage  of  all  structure  types  which  pass 
through  a  given  plane,  such  that  the  sum  of  all  of  the  structure  pdfs  produces  the  trajectory  velocity  pdf  for 
that  plane. 

In  Fig.  50a,  a  majority  of  the  model  structure  detected  are  attached.  Shorter  SCS  (i.e. ,  4/(1, 2))  convect 
slower  than  the  local  mean  (dashed  line  with  +  symbol)  on  average,  while  taller  SCS  (i.e.  4/(1, 4))  constitute 
the  part  of  the  pdf  traveling  faster  than  the  local  mean,  on  average.  For  each  attached  structure  type,  the 
ratio  of  the  mean  trajectory  velocity  to  the  local  mean  velocity  in  the  uppermost  plane  through  which  the 
structure  passes,  j,  ranges  from  0.85-0.88.  This  suggests  that  attached  hairpin  vortex  structures  identified 
by  SCS  do  not  convect  at  the  speed  of  their  heads,  on  average,  but  at  the  local  mean  somewhere  between 
their  head  and  the  wall. 

Further  from  the  wall,  as  in  Fig.  50b,  even  the  tallest  attached  structures  travel  at  a  speed  around  85%  of 
the  free-stream  velocity  ( u+  =  18.7)  on  average,  and  thus  the  high  speed  end  of  the  pdf  is  representative  of 
detached  structures.  Detached  structures  which  reach  down  closer  to  the  wall,  4/(2, 4)  and  4/(3, 4),  primarily 
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contribute  to  the  low-speed  end  of  the  pdf  while  structures  reaching  further  from  the  wall,  ^(4,4)  primarily 
contribute  to  the  high-speed  end. 


5.7  Summary  of  TRPIV  results 

We  presented  for  the  first  time  3D  streamwise  velocity  power  spectra  (2D  in  space  and  ID  in  time)  in  all 
regions  of  the  turbulent  boundary  layer  and  outlined  the  method  for  calculating  and  normalizing  these  from 
PIV  measurements.  The  effect  of  the  PIV  interrogation  window  on  the  resolution  and  attenuation  of  both  the 
spectrum  and  the  measured  velocity  fluctuations  was  discussed.  In  addition,  the  methods  used  for  reducing 
aliasing  and  spectral  leakage  were  outlined.  From  the  3D  spectra  presented,  the  more  usual  ID  and  2D 
spectra  were  also  calculated  and  compared  to  previous  findings.  The  3D  spectra  present  input  conditions 
for  the  resolvent  analysis  to  be  presented  subsequently. 

With  the  current  data  set,  it  was  found  that  the  spanwise  resolution  of  the  streamwise  velocity  fluctu¬ 
ations  was  lacking,  while  the  streamwise  and  temporal  resolution  were  reasonable.  Current  experiments 
are  underway  to  improve  the  spanwise  resolution  and  also  to  increase  the  dynamic  range  to  improve  the 
resolution  of  spanwise  fluctuations  so  that  the  spanwise  and  streamwise-spanwise  cross  spectrum  can  be 
analyzed.  Regardless  of  these  issues,  it  was  still  possible  to  consider  the  convection  velocity  of  streamwise 
scales  over  both  streamwise  and  spanwise  wavenumbers  using  the  calculated  3D  spectrum.  The  investigation 
and  interpretation  of  these  is  an  area  of  current  interest. 

Swirling  coherent  structures  (SCS)  were  reliably  extracted  from  two-dimensional  information  in  wall  parallel 
planes  and  tracked  in  time  with  geometric  properties  matching  other  data  reported  in  the  literature.  Using 
time-resolved  measurements,  several  important  conclusions  were  drawn.  First,  differences  were  found  in 
the  distribution  of  trajectory  velocities,  <  Utraj  >>  compared  to  the  average  core  velocity  calculated  from 
instantaneous  snapshots,  <  Ucore  >,  a  measure  often  reported  in  the  literature.  The  major  difference  was  a 
longer  low-speed  tail  using  the  trajectory  velocity,  indicating  that  the  time-resolved  data  and  the  centroid¬ 
tracking  algorithm  allowed  the  detection  of  slow  traveling  vortices  (although  the  error  in  <  Utraj  >  associated 
with  an  unquantified  potential  distribution  of  wall-normal  convection  velocities  is  noted).  Besides  the  low 
speed  tail,  the  distributions  from  both  calculations  were  nearly  identical. 

Time-resolved  measurements  also  allowed  the  measurement  of  the  lifespan  of  vortices,  which  increased  with 
increasing  wall-normal  distance,  and  a  small  percentage  of  vortices  persisted  for  more  than  5  Te  at  y/5  =  0.48. 
This  is  the  first  time  such  a  long  tracking  history  has  been  presented  in  a  wall-bounded  shear  flow.  This  is 
not  as  long  as  the  structure  lifetime  predicted  by  Elsinga  and  Marusic  [2010],  but  it  should  be  noted  that 
this  measure  was  inferred  from  a  study  of  the  rate  of  change  of  the  velocity  gradient  invariants,  Q  and  R, 
which  were  integrated  to  recover  trajectories  in  the  (Q,  R)  plane.  These  trajectories  were  found  to  form 
an  orbit  around  the  origin,  oscillating  between  the  4  basic  structure  types  (stable  and  unstable  nodes  and 
vortices),  where  the  time  to  complete  a  full  orbit  was  taken  as  a  measure  of  the  average  lifetime  of  coherent 
structures.  Considering  only  the  time  where  the  flow  topology  around  a  fluid  particle  would  be  a  stable 
focus,  reported  to  be  39%  of  the  orbit  time  in  a  turbulent  boundary  layer  by  Elsinga  and  Marusic  [2010], 
the  estimated  lifespan,  5.58  Te ,  would  be  closer  to  the  present  measurements.  It  should  also  be  reiterated 
that  since  the  current  measurements  were  taken  in  a  2D  plane,  it  is  possible  for  vortex  structures  to  have 
wall-normal  convection  velocities,  implying  that  the  measured  SCS  may  move  out  of  the  plane  of  view  and 
leading  to  lifetimes  that  may  in  some  cases  be  shorter  than  that  of  the  real  vortex  structures,  and  that  the 
streamwise  convection  velocity  pdfs  may  be  erroneously  broadened. 

Further  analyzing  the  time-resolved  data,  SCS  trajectories  were  found  to  be  streamwise  aligned  for  the 
most  part,  with  trajectory  angles  no  more  than  ±10°,  decreasing  with  increasing  distance  from  the  wall. 
This  finding  along  with  the  observations  of  the  evolution  of  the  full  time-resolved  velocity  fields  supports 
the  hypothesis  that  vortex  packets  do  not  grow  continuously  in  the  spanwise  direction  as  they  convect 
downstream,  but  must  do  so  in  some  discrete  manner. 

Finally,  using  the  convection  velocity  histograms,  and  under  the  limiting  assumption  that  a  coherent  structure 
in  the  flow  must  travel  at  or  near  a  single  velocity  at  all  distances  from  the  wall,  the  distribution  of  “model 
structures”  were  inferred  with  a  focus  on  whether  or  not  these  structures  were  attached,  that  is,  have 
a  signature  in  the  plane  nearest  the  wall.  From  joint  convection  velocity,  lifetime  pdfs,  it  was  found  that 
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attached  structures  decay  more  than  twice  as  fast  as  detached  structures  in  time.  In  addition,  the  convection 
velocity  pdf  for  each  structure  type  was  used  to  highlight  the  composition  of  the  convection  velocity  pdf 
with  SCS  traveling  both  faster  and  slower  than  the  local  mean.  Near  the  wall,  most  detected  structures 
were  attached  to  the  wall  and  there  was  a  large  spread  in  convection  velocity  depending  on  the  wall-normal 
extent  of  the  structure  in  question.  Further  from  the  wall,  the  pdf  was  dominated  by  detached  structures 
and  the  high  speed  tail  was  thinned  compared  to  the  low  speed  tail  because  structures  do  not  appear  to 
travel  faster  than  the  free-stream  velocity.  While  this  method  provided  results  consistent  with  the  attached 
eddy  hypothesis  and  capable  of  explaining  at  least  some  of  the  characteristics  of  SCS,  further  examination 
of  the  three-component,  three-dimensional  velocity  field,  either  from  experiments  or  simulations,  would  add 
confidence  in  the  importance  of  the  findings  described  here;  however  this  work  underlines  the  importance  of 
time-resolved  measurements  in  determining  the  characteristics  of  coherent  structure. 
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6  2D/3C  model  for  the  mean  velocity  in  turbulent  pipe  flow 


The  stochastically-forced,  two-dimensional,  three- (velocity)  component  model  for  the  mean  velocity  in  wall 
turbulence  was  introduced  by  Gayme  et  al.  [2010],  who  considered  channel  flow  in  earlier  work.  Under  the 
PECASE  award,  this  model  was  extended  to  pipe  flow.  The  results  revealed  some  important  differences 
between  pipe  and  channel  flow,  in  that  the  present  work  identified  the  footprint  of  the  turbulent  puffs  that 
are  observed  in  pipe  flow  transition.  The  capture  of  the  blunting  of  the  laminar  mean  velocity  towards  the 
turbulent  profile  is  an  important  first  step  to  providing  a  mean  profile  for  input  to  the  resolvent  analysis 
without  resorting  to  experimental  or  computational  results. 

The  work  described  below  was  published  as  Bourguignon  and  McKeon  [2011]  in  Physics  of  Fluids. 

6.1  Model  and  numerical  methods 

The  streamwise-constant  model  of  turbulent  pipe  flow  is  derived  from  the  NSE  written  in  cylindrical  coor¬ 
dinates  under  the  assumption  of  streamwise  invariance,  i.e. ,  it  constitutes  a  projection  of  the  NSE  onto  the 
streamwise  direction.  We  employ  a  nondimensionalization  based  on  the  pipe  radius  R  and  the  bulk  velocity 
U,  i.e.,  =  r/R ,  r  =  Ut/R ,  and  Re  =  2 RU/v,  where  v  is  the  kinematic  viscosity,  throughout  this  paper. 

Continuity  is  enforced  via  the  introduction  of  a  dimensionless  streamfunction  W  whose  evolution  equation 
is  obtained  by  taking  the  curl  of  the  NSE  projected  in  the  axial  direction.  The  model  consists  of  a  forced 
evolution  equation  for  the  streamfunction,  from  which  the  radial  and  azimuthal  velocities  can  be  derived,  and 
an  evolution  equation  for  the  axial  velocity  (in  terms  of  the  deviation  from  the  laminar  profile)  corresponding 
to  the  axial  momentum  balance,  subject  to  boundary  conditions  (BCs)  of  no-slip  and  no-penetration  on  the 
wall  of  the  pipe.  The  deviation  of  the  local  axial  velocity  from  laminar  illustrates  how  the  flow  evolves  away 
from  the  laminar  state  and  is  defined  as 


ux  =  ux(r})-U(ri),  (32) 

where  ux  is  the  instantaneous  axial  velocity  and  U ( r ])  =  1  -  r]2  is  the  laminar  base  flow. 

The  2D/3C  model  was  first  derived  by  Joseph  and  Tao  (1963)6  and  is  written  as  follows  for  the  cylindrical 
coordinate  system  shown  in  Figure  1: 
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Only  the  streamfunction  equation  is  forced,  based  on  the  results  of  the  study  by  Jovanovic  and  Barnieh 
[2005]  which  showed  that  maximum  amplification  is  obtained  by  forcing  in  the  cross-sectional  plane  in  the 
linearized  NSE.  Thus,  Afw  represents  a  forcing  term  that  is  required  to  maintain  the  perturbation  energy  in 
an  otherwise  stable  system  and  can  be  considered  to  represent  “noise  that  is  always  present  in  experiments, 
e.g.,  wall  roughness,  vibrations,  non-alignment  of  the  different  sections  of  the  pipe,  thermal  effects,  as  well  as 
taking  into  account  the  effects  not  modeled  by  the  streamwise  invariance  approximation.  In  the  subsequent 
sections,  we  consider  two  of  the  simplest  possible  forms  for  J\TW  in  order  to  investigate  the  origin  of  the 
blunting  of  the  mean  velocity  profile.  The  nonlinear  terms  in  the  governing  equation  for  the  streamfunction 
are  neglected  in  order  to  obtain  the  simplest  model  able  to  capture  the  blunting  of  the  velocity  profile  and 
also  because  their  effects  can  be  incorporated  into  the  unstructured  forcing  term  J\fv, .  Moreover,  the  study 
of  Gayme  et  al.  [2010]  showed  that  there  are  no  significant  differences  in  the  Couette  flow  statistics  obtained 
from  the  model  based  on  a  linearized  stream-  function  equation  compared  to  the  fully  nonlinear  2D/3C 
model.  The  bulk  velocity  is  maintained  constant  by  adjusting  the  pressure  gradient  C,  i.e.,  the  Reynolds 
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Figure  51:  The  coordinate  system  used  to  project  the  Navier-Stokes  equations. 


number  is  held  constant  for  each  study.  The  streamwise  velocity  behaves  as  a  passive  scalar  convected  by 
the  in-plane  velocities. 

The  2D/3C  model  with  stochastic  forcing  is  discretized  using  a  spectral-collocation  method  based  on  Cheby- 
shev  polynomials  in  the  radial  direction  and  Fourier  modes  in  the  azimuthal  direction,  associated  with  a 
third-order  semi-implicit  time  stepping  scheme  described  in  Spalart  et  al.  [1991].  The  singularity  at  the  origin 
of  the  polar  coordinate  system  is  avoided  by  re-defining  the  radius  from  -  1  to  1  and  using  an  even  number  of 
grid  points  in  the  radial  direction  [Heinrichs,  2004].  Three  Sylvester  equations  are  written,  respectively,  for 
A'l',  4'  and  ux,  associated  with  homogeneous  Dirichlet  boundary  conditions  and  are  solved  using  a  Fortran 
code  relying  on  an  optimized  Sylvester  equation  solver  from  the  SLICOT  numerical  library  [Jonsson  and 
Kagstrom,  2003].  The  BCs  4/  =  0  and  ux  =  0  at  the  wall  correspond,  respectively,  to  no-penetration  and 
no-slip  in  the  axial  direction.  The  no-slip  BC  in  the  azimuthal  direction  is  enforced  by  adding  particular 
solutions  to  the  streamfunction,  follow-  ing  the  influence  matrix  method  for  linear  equations  [Peyret,  2002], 
such  that  the  azimuthal  velocity  u $  =  -d'k/dr]  vanishes  at  the  wall.  Under  deterministic  forcing,  the  2D/3C 
model  is  reduced  to  a  set  of  two  ordinary  differential  equations  (ODEs)  that  are  solved  in  MATLAB  using 
spectral  methods  based  on  a  Chebyshev  polynomial  expansion. 


6.2  Simplified  2D/3C  model  with  deterministic  forcing 


We  begin  by  developing  a  simplified  version  of  the  2D/3C  model  subject  to  a  steady,  deterministic  forcing 
to  study  momentum  transfer  between  the  in-plane  and  axial  velocities.  The  study  of  optimal  disturbance 
growth  in  pipe  flow  by  Schmid  and  Henningson  [2001]  demonstrated  that  the  streamwise-constant  mode  with 
azimuthal  wavenumber  n  =  1  is  the  most  amplified  based  on  an  energy  norm.  Thus,  we  isolate  this  mode 
as  a  candidate  perturbation  contributing  to  the  blunting  of  the  velocity  profile  and  consider  a  forcing  with 
only  this  one  mode  in  the  azimuthal  direction,  namely 


=  N(if)  sin^.  (36) 

The  streamfunction  SP (77, 0)  has  the  same  azimuthal  dependence  as  the  forcing  since  its  governing  equation 
is  linear  i.e.,  4'  =  B/ 1  (77)  sin^.  The  axial  velocity  can  be  written  in  terms  of  a  mean  deviation  from  laminar 
uq  and  a  zero-mean  perturbation  u±  cos  </>  corresponding  to  the  linear  response  of  the  system  to  the  forcing, 

ux{ri,(j>)  =  ^0(77)  +Mi(?y)cos^.  (37) 


The  2D/3C  model  can  be  simplified  as  follows  to  predict  the  steady  state  mean  deviation  from  laminar,  uq, 
obtained  with  the  deterministic  forcing  profile  N(rj): 


(A,—)  ^  =  -0.5 ReN{r,), 
^A,;  —  j  ui  =  0.5  Re'S?  idv(U  +  uq), 


(38) 

(39) 
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A,;mo  = -0.5i?e(C?7  -  (40) 

where  =  r/Dvv  +  dv  is  the  radial  derivative  component  of  the  2D  Laplacian  premultiplied  by  77  and 
Dv  =  d / d?7 .  In  order  to  obtain  the  simplest  model  able  to  capture  the  blunting  of  the  velocity  profile,  we  can 
linearize  equation  40  under  the  assumption  of  small  amplitude  forcing.  The  resulting  model  contains  only 
one  nonlinear  term  in  one  ODE,  the  other  two  ODEs  being  linear.  The  presence  of  at  least  one  nonlinear 
term  is  required  to  obtain  a  change  in  mean  flow  since  linear  models  always  give  the  same  mean  flow  as  the 
one  used  for  the  linearization. 

A  simple  inspection  of  equations  38-40  leads  to  the  following  observations.  Conclusions  similar  to  those  of 
Reddy  and  Ioannou  [2000]  on  the  energy  transfer  between  streamwise  constant  modes  can  be  recovered,  but 
this  time  in  terms  of  momentum  transfer  and  for  the  pipe  instead  of  Couette  flow.  The  mean  turbulent 
mode  uo  cannot  extract  momentum  from  the  laminar  base  flow  and  is  sustained  by  the  nonlinear  interaction 
between  the  axial  velocity  perturbation  u\  and  the  streamfunction  Ti,  i.e. ,  the  term  in  equation  40. 

The  non-normality  of  the  system  manifests  itself  in  the  linear  coupling  between  the  laminar  base  flow  and 
the  streamfunction  which  amplifies  the  disturbances  and  generates  the  axial  velocity  perturbation  iq  by 
convection  of  streamwise  momentum  (see  equation  39).  The  shape  of  the  streamfunction  determines  the 
amount  of  blunting  obtained  for  a  given  amplitude  coefficient. 

In  order  to  advance  further  analytically,  we  write  the  streamfunction  profile  f !  as  a  Taylor  series  at  the 
origin,  i.e., 

=  ££o<W,  (41) 

and  set  ao  =  0  in  order  to  enforce  continuity  in  the  limit  of  77  tending  to  zero,  recalling  that  (77, 0)  = 
-ty{j],(j)  +  7 r).  The  forcing  profile  generating  Ti  is  given  by 

=  (42) 

We  rescale  the  coefficients  cq  by  aq,  i.e.,  =  aq[?7  +  a2?72  +  CK3773  +  aq 77 4  +  ...]  and  choose  oq  such  that  the 

change  in  mean  flow  induced  by  the  forcing  has  the  same  amplitude  at  its  maximum  as  in  the  experiments 
of  den  Toonder  and  Nieuwstadt  [1997]  at  the  same  Reynolds  number,  in  order  to  facilitate  comparison  of 
the  results.  Note  that,  depending  on  the  streamfunction  profile,  the  coefficient  al  is  not  necessarily  small 
in  which  case  the  linearization  of  equation39  is  no  longer  justified.  In  the  following,  we  solve  the  nonlinear 
momentum  balance  for  u\  (Eq.  (6))  regardless  of  the  amplitude  of  oq.  A  fundamental  streamfunction  profile 
Ti  =  aq(?7  -  3?73  +  2?74)  is  obtained  by  truncating  the  series  expansion  to  the  fourth  order  term,  enforcing  the 
BCs  4/i  =  city  1 /dr/  =  0  at  the  wall,  and  requiring  that  the  forcing  be  bounded  at  the  origin. 

Streamfunctions  given  by  Tl  a  =  Tj,  Ti  t,  =  Tf  and  Ti  c  =  were  investigated  in  order  to  ascertain  the 
ability  of  such  simple  functions  to  capture  key  aspects  of  the  blunting  of  the  mean  profile  and  to  identify 
the  role  of  the  radial  streamfunction  profile.  The  amplitude  coefficients,  oq,  were  chosen  such  that  the 
same  amount  of  blunting  is  realized  in  each  case,  as  described  above  in  terms  of  the  maximum  deviation 
from  laminar.  The  steady-state  equations  were  solved  with  64  grid  points  in  the  radial  direction,  and  at 
a  Reynolds  number  of  24  600,  matching  one  of  the  pipe  flow  experiments  of  den  Toonder  and  Nieuwstadt 
[1997].  A  short  convergence  study  showed  that  this  resolution  in  the  radial  direction  is  sufficient  since  the 
maximum  relative  error  com-  pared  to  the  solution  computed  on  192  grid  points  is  less  than  one  percent. 

Figure  2  shows  the  radial  forms  of  the  three  analytic  streamfunctions,  and  the  respective  resulting 

variations  of  the  mean  deviation  from  laminar.  Streamfunction  profiles  4/i ^  and  TiiC  reach  maximum 
amplitudes  of  about  0.05  and  0.25,  respectively,  in  the  core  of  the  pipe,  while  in  comparison  HT  a  is  relatively 
flat  with  a  maximum  amplitude  of  0.0085.  Despite  such  a  wide  variation  in  streamfunction  amplitude  between 
the  three  cases,  the  streamwise  velocity  profiles  are  remarkably  similar.  Even  the  simple  streamfunction 
profile  Tqa  leads  to  a  good  blunting  of  the  velocity  profile,  in  the  sense  that  the  general  features  of  the  mean 
pro-  file  are  reproduced.  The  maxima  of  the  velocity  profiles  are  situated  further  from  the  wall  compared 
to  the  experimental  data  [den  Toonder  and  Nieuwstadt,  1997]  which  likely  corresponds  to  the  neglect  of  the 
influence  of  the  small  scales  near  the  wall  by  the  streamwise-constant  model.  The  results  obtained  with  the 
simplified  2D/3C  model  also  show  that  the  velocity  profile  is  relatively  independent  of  the  radial  shape  of 
the  forcing  and  streamfunction,  i.e.,  the  profile  can  be  said  to  be  robust  to  the  shape  of  the  streamfunction. 
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Figure  52:  (a)Streamfunctions  ^ri>0-c(»7)  and  (b)  corresponding  velocity  profiles  uo(r])  for  4 ’i,a(v)  = 
0.033(77  -  3?73  +  2t74)  (thin  solid),  =  0.7(7/  -  3r f  +  2r /4)2  (dashed),  and  ^Efi,c(77)  =  14(ry  -  Sr]3  +  2r]4)3 

(dash-dot)  and  experimental  velocity  profile  of  den  Toonder  and  Nieuwstadt  [1997]  at  Re  =  24600  (thick 
solid). 


The  in-plane  kinetic  energy,  defined  as  the  integral  of  m2/2  +  'u2/2  over  the  pipe  cross  section,  varies  from 
3.7a;10-4  for  4/la  to  0.31  for  4T  c  even  though  the  same  amount  of  blunting  is  realized  by  the  three  stream- 
functions  4T  ac.  This  large  variation  of  the  in-plane  kinetic  energy  between  different  streamfunctions  (about 
3  orders  of  magnitude)  can  be  understood  by  consideration  of  the  influence  of  the  near-wall  region  in  the 
pipe  on  the  overall  amplification.  Amplification  is  proportional  to  the  mean  shear  and  surface  area,  which 
are  both  maximum  at  the  wall.  Hence,  the  streamfunction  47ic  having  a  small  amplitude  near  the  wall 
compared  to  see  Figure  52,  does  not  benefit  from  the  near-wall  effects  as  much  as  l'ii0  and  needs  to 
reach  a  larger  amplitude  in  the  core  (about  30  times  larger)  in  order  to  give  the  same  amount  of  blunting  as 
47 1  a,  resulting  in  a  significantly  larger  in-  plane  kinetic  energy. 

The  large  amplitude  reached  by  47i>c  in  the  core  of  the  pipe  results  in  a  maximum  radial  velocity  amplitude 
of  64%  of  the  bulk  velocity  compared  to  3.3%  for  Tiia,  significantly  larger  than  the  radial  velocity  turbulence 
intensity  measured  by  den  Toonder  and  Nieuwstadt  [1997]  which  is  about  1  plus  unit  or  equivalently  about 
5%  of  the  bulk  velocity.  A  streamfunction  that  maintains  a  forcing  amplitude  comparable  to  the  experimental 
noise  level  should,  therefore,  have  a  larger  amplitude  near  the  wall  and  relatively  constant  amplitude  over  the 
whole  domain.  A  similar  conclusion  can  be  obtained  by  considering  that  blunting  results  from  the  advection 
of  axial  momentum  by  the  radial  velocity  so  that  a  large  amount  of  blunting  is  realized  when  large  radial 
velocities  are  present.  Taking  into  account  that  the  radial  velocity  depends  on  the  azimuthal  velocity  via 
the  continuity  constraint,  it  can  be  seen  that  maximization  of  the  ratio 

suggests  that  the  flattest  streamfunction  profile  (or  equivalently  the  simplest  radial  dependence  results  in 
the  largest  amplification.  In  terms  of  structures,  the  largest  structures  corresponding  to  the  modes  with  the 
least  zero  crossings  in  the  radial  direction  are  more  able  to  redistribute  momentum  over  the  cross-section 
of  the  pipe.  The  impor-  tance  of  the  modes  with  the  least  zero  crossings  is  a  known  feature  of  turbulent 
pipe  flow:  modes  with  a  radial  quantum  number  of  1  in  the  study  of  Duggleby  et  al.  [2007]  were  shown  to 
capture  most  of  the  energy  in  their  dynamical  eigenfunction  decomposition  of  turbulent  pipe  flow.  Likewise, 
the  singular  modes  that  are  most  amplified  in  the  study  of  McKeon  and  Sharma  [2010]  exhibit  the  lowest 
number  of  zero  crossings. 

Based  on  our  numerical  study,  the  streamfunction  47(77)  =  0.033(?7  —  3rj3  +  2rj4)  sirup  leads  to  a  blunting  of  the 
velocity  profile  whose  maximum  amplitude  matches  the  experimental  data  at  the  same  Reynolds  number 
and  is  generated  by  a  forcing  profile  N(r 7)  =  -90/ Re  that  is  constant  in  the  radial  direction  and  whose 
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(a)  (b)  (c) 

Figure  53:  Model  output  for  deterministic  forcing:  (a)  contours  of  the  streamfunction  4/1>a(/y)  =  0.033(?y- 
3?y3  +  2 ry4)  sin  0,  (b)  vector  plot  of  the  corresponding  in-  plane  velocities,  and  (c)  contours  of  the  resulting 
axial  velocity  Held. 


amplitude  over  the  pipe  cross-section  is  consistent  with  the  experimentally  measured  rms  amplitude  of  the 
turbulence  fluctuations.  The  contours  of  the  streamfunction  T(?y)  =  0.033(ry  -  3 ry3  +  2 rf)  sin^  are  plotted  on 
Figure  53  together  with  a  vector  plot  of  the  corresponding  in-plane  velocities  and  the  resulting  contours  of 
the  axial  velocity  field.  The  streamfunction  exhibits  two  counter-rotating  rolls  which  advect  the  mean  shear 
to  create  a  low-  and  a  high-speed  streak  of  axial  velocity  defined  with  respect  to  the  laminar  base  flow.  The 
high-speed  streak  sits  near  the  wall  whereas  the  low-speed  streak  is  localized  near  the  centerline.  Hence, 
the  flow  is  on  average  faster  near  the  wall  and  slower  at  the  center  compared  to  laminar  as  is  the  case  for 
the  velocity  profile  of  turbulent  pipe  flow.  The  amplification  factor  between  the  in-plane  and  streamwise 
velocities,  defined  as  the  ratio  of  the  extrema,  is  about  20  in  this  case. 

As  a  final  comment  related  to  this  simple,  steady,  deterministic  streamfunction  analysis,  we  note  that  when 
the  azimuthal  wavenumber  is  chosen  to  match  the  azimuthal  dependence  of  the  spatial  puffs  observed  by  Hof 
et  al.  [2004]  and  a  new  streamfunction  model  is  derived  along  the  lines  described  above,  the  velocity  fields 
produced  by  our  model  are  remarkably  similar  to  the  measured  velocity  fields  near  the  trailing  edge  of  the 
puffs.  For  an  azimuthal  wavenumber  equal  to  six,  the  lowest  order  streamfunction  profile  that  satisfies  the 
BCs  and  is  generated  by  a  forcing  profile  bounded  at  the  origin  is  given  by  4/6  =  aq(?y4  ~  2ry5  +  ry6);  figure  4 
shows  the  resulting  variation  of  the  axial  velocity.  The  wall-normal  position  of  the  high-  and  low-speed  axial 
velocity  streaks  compares  well  with  figures  2(e)  and  2(f)  in  the  paper  by  Hof  et  al.  [2004].  Moreover,  our 
model  appears  to  capture  the  merging  of  the  low-speed  streaks  and  their  congregation  near  the  center  of  the 
pipe  which  is  observed  in  experiments  [Hof  et  al.,  2004]  but  is  not  present  in  the  traveling  wave  solutions  of 
the  NSE. 

Despite  the  simplicity  of  the  deterministic  streamfunction  profiles  described  in  this  section,  key  aspects  of 
the  axial  velocity  variation  observed  in  experiments  are  recovered.  In  Sec.  IV,  we  consider  a  more  realistic, 
time-dependent  forcing  function  in  an  effort  to  capture  more  details  associated  with  the  mechanisms  of 
momentum  transfer. 


6.3  Stochastic  forcing  of  the  2D/3C  model 

Effects  that  are  not  captured  by  the  (unforced)  2D/3C  model  and  external  perturbations  that  are  always 
present  in  experiments  but  not  explicitly  described  by  the  NSE  are  likely  much  better  captured  by  a  stochastic 
forcing  than  a  deterministic  one.  In  this  section,  we  present  the  results  from  a  simulation  of  the  2D/3C  model 
forced  by  small  amplitude  white  noise,  an  approach  which  was  successfully  explored  by  Gayme  et  al.  [2010] 
in  Couette  flow  as  well  as  several  previous  linear  studies.  Such  stochastic  forcing  has  the  advantage  over 
deterministic  forcing  of  not  relying  on  any  assumption  regarding  the  spatial  and  temporal  dependence  of  the 
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Figure  54:  Contours  of  the  axial  velocity  induced  by  the  streamfunction  <j>)  =  (rf  -  2r]5  +  rf)  sin(6^), 
the  light  and  dark  filled  contours  correspond  to  regions  of  the  flow,  respectively,  faster  and  slower  than 
laminar. 

perturbations,  such  that  the  resulting  velocity  field  reflects  the  direction  of  maximum  (assumed  dominant) 
disturbance  amplification.  The  noise,  iV^,,  is  applied  at  each  grid  point  in  space  and  at  every  time  step 
and  follows  a  normal  distribution  with  zero  mean  and  a  variance  that  depends  on  the  radius  such  that  the 
variance  per  surface  area  is  constant.  In  order  to  prevent  aliasing  in  the  nonlinear  coupling  terms  in  the 
streamwise  velocity  equation,  we  truncate  the  2D  Fourier  transform  of  the  forcing  term  after  the  lowest 
two-third  wavenumbers,  as  described  in  Canuto  [2006]. 

Representative  time  traces  of  the  centerline  velocity  for  two  different  Reynolds  numbers  are  reproduced  on 
Figure  5  and  show  numerous  sharp  drops,  which  we  identify  as  the  signature  of  streamwise-constant  puffs, 
before  increasing  smoothly  nearly  back  to  its  laminar  value.  (Since  there  is  no  grid  point  at  the  centerline, 
we  approximate  the  centerline  velocity  by  averaging  the  axial  velocity  in  the  azimuthal  direction  over  the 
grid  points  closest  to  the  center  of  the  pipe.)  The  signatures  of  the  streamwise-constant  puffs  are  remarkably 
similar  to  the  spatial  evolution  of  the  centerline  velocity  from  the  trailing  to  the  leading  edge  of  the  spatial 
puffs  in  the  numerical  simulations  of  Shimizu  and  Kida  [2009] .  We  define  a  puff  generation  timescale  as  the 
time  elapsed  between  two  sharp  drops  of  the  centerline  velocity.  At  Re  =  2200,  the  timescale  is  about  75 
dimensionless  time  units  based  on  the  pipe  radius  and  compares  well  with  the  time  scale  of  puffs  computed 
from  experimental  data[Nishi  et  al.,  20078],  in  which  the  puffs  are  5  to  20  diameters  long  depending  on 
Reynolds  number  and  convect  at  nearly  the  bulk  velocity.  A  puff  of  length  20D  which  is  separated  from  the 
next  puff  by  a  laminar  region  of  length  equivalent  to  one  puff  would  lead  to  a  dimensionless  timescale  based 
on  the  pipe  radius  of  80. 

The  puff  generation  timescale  is  an  increasing  function  of  the  Reynolds  number,  reflecting  the  fact  that 
puffs  tend  to  be  longer  in  experiments  as  the  Reynolds  number  increases  Figures  5(a)  and  5(b),  reaching 
330  at  Re  -  10000,  but  is  relatively  independent  of  the  noise  amplitude  Figures  5(a)  and  5(c).  The  drop  in 
centerline  velocity  associated  with  the  trailing  edge  of  a  puff  can  be  observed  to  be  sharper  and  stronger  for 
larger  forcing  amplitudes.  Note  that  if  we  relax  the  BCs  to  allow  for  slip  in  the  azimuthal  direction,  i.e. , 
we  use  a  shear-stress  free  condition  (not  shown) ,  the  simulations  capture  the  creation  of  streamwise  vortices 
and  streaks  as  well  as  the  blunting  of  the  velocity  profile,  but  we  do  not  observe  clearly  the  cyclic  generation 
of  puffs  in  the  time  evolution  of  the  full  velocity  field  or  their  signature  in  the  time  traces  of  the  centerline 
velocity.  However,  for  a  given  forcing  amplitude,  the  amount  of  blunting  realized  with  the  slip  BC  in  the 
azimuthal  direction  is  larger  than  with  no-slip. 

The  time  evolution  of  the  flow  field  is  characterized  by  the  quasi-periodic  generation  of  streamwise-constant 
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Figure  55:  Time  traces  of  the  centerline  velocity  from  three  different  simulations,  respectively,  at  Re  =  2200 
with  0.0005  and  0.002  rms  noise  levels  (a)  and  (c)  and  at  Re  =  10000  with  0.002  rms  noise  level  (b).  The 
resolution  in  the  radial  direction  is  N  =  48.  (d)  Zoom  on  the  time  interval  during  which  the  samples  of 
Figure  6  are  taken.  The  vertical  lines  indicate  the  sampling  instants. 


puffs  followed  by  their  decay  and  the  return  of  the  flow  close  to  the  laminar  state,  i.e. ,  each  bursting  event 
is  followed  by  quiescent  flow  equivalent  to  the  laminar  regions  that  separate  the  puffs  in  the  experiments 
of  Lindgren  [1957].  The  three  main  stages  in  the  evolution  of  the  flowheld  corresponding  to  a  streamwise- 
constant  puff  are  plotted  on  Figure  6  in  terms  of  the  axial  velocity  ((a)  to  (c))  and  the  swirling  strength 
((d)  to  (f)),  defined  as  the  magnitude  of  the  imaginary  part  of  the  in-plane  velocity  gradient  eigenvalues 
and  representative  of  coherent  vorticity  without  the  influence  of  mean  shear.  These  instantaneous  fields 
correspond  to  the  time  instants  marked  with  vertical  lines  in  Figure  5(d).  During  the  first  stage,  patches 
of  swirl  move  toward  the  center  of  the  pipe  Figures  6(d)  and  5(e)  and  create  streaks  by  convection  of  the 
axial  momentum  Figure  6(a).  The  radial  motion  of  the  coherent  swirl  corresponds  to  a  lift  up  of  streamwise 
vortices  away  from  the  wall  as  observed  in  the  simulations  of  van  Doorne  and  Westerweel  [2009]  if  we  consider 
the  evolution  of  the  vorticity  field  projected  in  a  plane  moving  at  the  bulk  velocity.  The  second  stage  consists 
of  the  segregation  of  the  high-  and  low-speed  streaks,  the  latter  being  convected  toward  the  center  of  the 
pipe  resulting  in  the  blunting  of  the  velocity  profile  characteristic  of  turbulent  pipe  flow,  Figure  6(b).  Once 
a  low-speed  streak  reaches  the  center,  the  centerline  velocity  drops  sharply,  as  can  be  seen  on  the  time  traces 
on  Figure  5(b).  Finally,  the  swirling  strength  and  streaks  decay  Figures  6(c)  and  6(f)  and  the  flow  returns 
close  to  the  laminar  state  before  the  next  cycle  starts.  As  the  streamwise  vortices  convect  toward  the  center 
of  the  pipe,  the  sign  of  the  net  azimuthal  velocity  is  reversed  such  that  the  rotation  sign  changes  from  one 
cycle  to  the  next. 

While  the  output  of  the  2D/3C  model  is  a  temporal  variation  of  a  streamwise  constant  field,  a  simplistic 
comparison  with  spatial  experimental  results  may  be  made  by  assuming  an  appropriate  convection  velocity 
to  be  the  bulk  velocity.  With  this  in  mind,  the  time  evolution  of  the  velocity  field  is  remarkably  similar  to 
the  flow  visualizations  by  Hof  et  al.  [2004]  in  transitioning  pipe  flow  when  a  puff  is  observed  in  a  reference 
frame  moving  with  the  bulk  velocity.  Streak  merging  in  experiments  was  reported  by  those  authors,  who 
showed  that  the  number  of  streaks  in  the  cross-section  decreases  due  to  their  merging  as  the  cross-stream 
observation  plane  is  moved  from  the  trailing  edge  to  the  leading  edge  of  a  puff.  This  streak  merging,  as  well 
as  the  segregation  of  the  high-  and  low-speed  streaks  observed  in  the  experiments  is  accurately  reproduced 
by  the  model,  as  shown  in  Figure  5.  The  time  evolution  of  the  vorticity  field  in  our  simulations  is  also 
consistent  with  the  projection  of  the  3-dimensional  vorticity  field  inside  the  puff  computed  by  van  Doorne 
and  Westerweel  [2009]  in  a  plane  moving  at  the  bulk  velocity,  characterized  by  a  lift  up  of  streamwise  vortices 
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(f) 


Figure  56:  Contours  of  the  axial  velocity,  subfigures  (a)  to  (c),  and  of  the  swirling  strength  for  the  in-plane 
velocities,  subfigures  (d)  to  (f),  computed,  respectively,  at  =  1620,  t  =  1700,  and  t  =  1740  dimensionless  time 
units. 


away  from  the  wall  as  we  move  from  the  trailing  to  the  leading  edge  of  the  puff  (corresponding  to  increasing 
time  in  our  model). 

The  generation  and  subsequent  decay  of  the  puffs  can  be  described  in  terms  of  a  simple  process  driven  by 
background  noise,  sketched  in  the  inner  loop  of  Figure  7.  In  essence,  the  presence  of  background  noise  in  the 
pipe  cross-section  results  in  the  formation  of  streamwise-constant  vortices  which  advect  axial  momentum  to 
create  high-  and  low-speed  streaks  of  axial  velocity.  The  vortices  interact  nonlinearly  with  the  streaks  to 
segregate  them,  i.e.,  to  convect  the  low-speed  streaks  toward  the  center  and  the  high-speed  streaks  toward  the 
wall,  leading  to  the  blunting  of  the  velocity  profile.  We  term  this  a  quasi-self-sustaining  process  (QSSP)  to 
reflect  the  fact  that,  while  there  is  no  feedback  from  the  axial  velocity  to  the  evolution  of  the  streamfunction 
(denoted  by  the  dashed  line  in  the  inner  loop  of  Figure  7)  which  could  sustain  the  cycle  described  in  Figure  6, 
the  insertion  of  small-amplitude  stochastic  forcing  in  the  cross-stream  plane  appears  to  provide  an  effective 
replacement  for  this  mechanism.  The  radial  shape  of  the  forcing  is  directly  related  to  the  amount  of  swirl 
present  in  the  simulations  and,  therefore,  suggests  a  method  to  control  the  flow  by  shaping  the  noise  forcing 
in  the  spirit  of  the  simple  control  mechanism  developed  by  Hof  et  al.  [2010]  which  reduces  the  inflection 
points  of  the  velocity  profile  and  leads  to  a  relaminarization  of  the  flow. 

A  well-studied,  fully  self-sustaining  process  has  been  described  for  the  (3D/3C)  NSE  by  Waleffe  [1997],  which 
relies  on  the  nonlinear  interaction  of  the  disturbances  developing  from  the  instability  of  the  streaks  to  force 
the  stream-  wise  vortices.  Our  results  suggest  that  the  blunting  of  the  velocity  profile  can  be  considered  to  be, 
in  essence,  a  nonlinear,  two-dimensional  phenomenon  in  which  the  directional  amplification  of  the  2D/3C 
model  selects  the  appropriate  disturbance  shape  from  stochastic  forcing  in  order  to  maintain  the  QSSP. 
The  model  suggests  that  the  sustaining  mechanisms  are,  therefore,  insensitive  to  the  type  of  (nonlinear) 
interaction  invoked  to  force  the  streamwise  vortices.  Based  on  the  simulations  presented  here,  we  argue  that 
it  is  the  presence  of  appropriate,  small  disturbances  in  the  flow  (necessarily  contained  in  the  white  noise 
of  stochastic  forcing),  and  not  their  specific  interaction  that  sustains  the  streamwise  vortices  and  streaks. 
Hence,  the  2D/3C  model  shows  that  using  unstructured  (stochastic)  forcing  to  generate  the  streamwise 
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at  the  wall  velocity 

Figure  57:  Diagram  detailing  the  different  stages  of  the  QSSP.  The  dashed  lines  represent  unmodeled 
effects. 

vortices  appears  to  be  sufficient  to  capture  the  formation  of  the  streaks  and  their  segregation  resulting  in  a 
blunting  of  the  velocity  profile  and  that  the  overall  flow  dynamics  are  relatively  insensitive  to  the  particular 
regeneration  process  invoked  to  produce  the  stream-wise  vortices.  The  QSSP  is  simpler  than  other  processes 
described  in  the  literature  but  does  appear  to  capture  the  minimum  turbulence  dynamics  and  produce  flow 
fields  dominated  by  streamwise  vortices  and  streaks  whose  temporal  evolution  and  topology  compares  well  to 
experimental  visualizations  and  numerical  simulations  of  puffs,  under  an  appropriate  projection  onto  a  plane 
convecting  at  the  bulk  velocity.  In  addition,  our  simulations  are  significantly  less  computationally  intensive 
than  studies  of  the  full  NSE[van  Doorne  and  Westerweel,  2009,  Shimizu  and  Kida,  2009]  since  our  domain 
is  2D  and  we  do  not  need  to  track  the  position  of  the  puffs  in  time. 

Interestingly,  our  simulations  of  the  stochastically  forced  2D/3C  model  for  pipe  flow  do  not  reach  fully 
developed  turbulence,  regardless  of  the  amount  of  forcing  and  Reynolds  number,  even  though  the  same 
model  applied  to  Couette  flow  and  described  in  Ref.  8  reached  a  steady  state  with  a  velocity  profile  in 
good  agreement  with  the  profiles  from  full  3D  simulations  at  the  same  Reynolds  number,  provided  that  the 
amount  of  forcing  is  appropriately  chosen.  Mellibovsky  et  al.  [2009]  describe  pipe  flow  transition  as  a  two- 
stage  process,  the  first  stage  consists  in  the  formation  of  the  puffs  and  the  second  in  their  spreading  in  space. 
Those  authors  argue  that  the  two  stages  are  caused  by  different  instability  mechanisms.  Our  simulations 
support  this  point  of  view  since  we  are  able  to  capture  the  first  stage  using  a  streamwise-constant  model  but 
not  the  second  stage  for  which  three-dimensional  effects  are  necessary  to  allow  for  the  puffs  to  destabilize 
and  form  slugs  which  will  expand  to  cover  the  whole  flow  domain. 

6.4  Summary 

We  have  investigated  a  (globally  stable)  streamwise-constant,  2D/3C  model  of  turbulent  pipe  flow  subject 
to  two  types  of  very  simple  forcing  profiles,  namely  a  deterministic,  steady  streamfunction,  and  stochastic 
forcing  of  the  stream-function  equation.  By  retaining  (only)  the  nonlinearity  coupling  the  (cross-plane) 
streamfunction  and  the  axial  velocity,  the  model  permits  the  study  of  the  physics  underlying  mo-  rnentum 
transfer.  This  model,  which  is  significantly  more  tractable  than  the  full  Navier-Stokes  equations,  appears 
to  be  capable  of  capturing  key  features  associated  with  transition  to  turbulence,  i.e.,  the  blunting  of  the 
velocity  profile  and  the  generation  of  streamwise  vortices  and  streaks. 

Using  the  simplest,  time- invariant  deterministic  forcing,  we  showed  that  the  velocity  profile  is  robust  with 
respect  to  variations  in  the  forcing  profile  and  we  produced  realistic  velocity  fields  that  are  remarkably  similar 
to  the  flow  visualizations  by  Hof  et  al.  [2004]  near  the  trailing  edge  of  a  puff.  Thanks  to  the  significant 
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reduction  in  complexity  associated  with  the  2D/3C  projection  compared  to  the  full  Navier-Stokes  equations, 
several  observations  can  be  made  by  examination  of  the  governing  equations:  momentum  is  extracted  from 
the  laminar  base  flow  by  (zero-mean)  perturbations  via  linear  non-normal  mechanisms  and  redistributed  via 
nonlinear  interactions  of  the  perturbations  to  result  in  a  change  in  mean  flow.  Hence,  our  simplistic  model 
allowed  us  to  isolate  basic  mechanisms  leading  to  the  blunting  of  the  velocity  profile  in  pipe  flow. 

Under  stochastic  forcing,  the  model  generates  “streamwise-constant  puffs”  at  a  frequency  that  depends  on 
the  Reynolds  number  but  not  on  the  forcing  amplitude.  The  model  captures  the  first  stage  of  pipe  flow 
transition  as  described  by  Mcllibovsky  et  al.  [2009],  i.e. ,  the  formation  of  puffs  but  not  the  second  stage 
(spreading  of  the  puffs  in  space)  for  which  some  three-dimensional  effects  are  necessary.  The  time  evolution 
of  the  velocity  fields  produced  by  our  simulations  is  remarkably  similar  to  flow  visualizations  in  transitioning 
pipe  flow  [Hof  et  ah,  2004]  when  a  puff  is  observed  in  a  reference  frame  moving  at  the  bulk  velocity.  The 
segregation  of  the  high-  and  low-speed  streaks  observed  in  experiments  is  accurately  captured  by  the  model 
as  well  as  the  streak  merging  and  the  lift  up  of  the  streamwise  vortices  away  from  the  wall. 

The  essential  dynamics  governing  the  generation  of  puffs  in  pipe  flow  transition  were  captured  by  the 
streamwise-constant  model  and  are  relatively  insensitive  to  the  particular  regeneration  mechanisms  invoked 
to  produce  the  streamwise  vortices,  permitting  the  introduction  of  a  quasi  self-sustaining  process  to  describe 
the  generation  of  puffs  and  the  blunting  of  the  mean  velocity  profile. 
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7  Conclusions  and  Outlook 


This  PECASE  award  has  enabled  the  development  of  resolvent  analysis  for  wall  turbulence,  a  promising 
approach  to  understanding,  and  ultimately  controlling  flows  of  interest  to  a  range  of  Air  Force  vehicles.  In 
addition,  experiments  and  modeling  have  been  performed  to  support  the  analysis  and  provide  appropriate 
input  parameters  and  profiles.  We  have  described  here  how  a  straightforward  analysis  permits  significant 
simplifications  in  complexity,  dimension,  and  interpretation  of  the  NSE.  The  most  important  of  these  may  be 
the  low  dimensionality  and  sparse  nature  of  the  formulation,  which  open  the  door  to  the  use  of  state-of-the- 
art  mathematical  tools  including  compressed  sampling  and  techniques  for  low  rank  matrix  approximation  of 
linear  operators  [Sharma,  2009].  A  study  in  a  turbulent  channel  [Moarref  et  ah,  2013],  enabled  by  an  efficient 
randomized  algorithm  for  the  approximation  of  the  SVD  of  low  rank  operators  [Halko  et  ah,  2011],  gives 
encouraging  results  with  respect  to  utilizing  the  low  rank  nature  of  the  resolvent  for  both  prediction  and 
modeling  purposes:  it  appears  that  a  very  good  approximation  to  the  distribution  of  streamwise  turbulent 
intensity  (the  dominant  contributor  to  the  turbulent  kinetic  energy)  can  be  obtained  using  only  the  first  two 
singular  modes  and  extraordinarily  simple  assumptions  about  the  amplitude  distribution  in  K  space.  The 
sparse  nature  of  the  modes  identified  as  being  active  in  real  flows  in  the  frequency  domainBourguignon  et  al. 
[2013]  naturally  leads  to  consideration  of  compressive  sampling  as  a  tool  for  data  reduction  [Bourguignon 
et  ah,  2014a].  Additionally,  the  resolvent  analysis  may  have  implications  for  simulation,  since  the  frequency 
domain  representation  effectively  captures  aspects  of  the  dynamics  that  are  masked  in  the  time  domain, 
i.e.  by  time-stepping  approaches.  This  grant  has  seeded  several  ongoing  collaborative  effort  aimed  at 
utilizing  resolvent  analysis  towards  control  of  wall  turbulence.  We  propose  that  the  success  of  the  resolvent 
analysis  described  here  for  incompressible,  Newtonian  flow  in  canonical  flow  configurations  suggests  promise 
for  generalization  to  more  complex  flows.  In  particular,  work  is  ongoing  to  consider  adaptation  of  the 
formulation  to  consider  rough-wall,  non-Newtonian  and  compressible  flows,  and  spatial  inhonrogeneity.  This 
will  facilitate  understanding  and  control  of  a  broad  range  of  practical  and  applied  flows. 
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